PYRENAIC ROCK GLACIERS: AN AIRBORNE AND MULTITEMPORAL LiDAR MONITORING CASE STUDY IN THE BESIBERRI AREA

The monitoring of rock glaciers is a current subject of interest because of its application as permafrost indicator and its sensitivity to climatic changes (especially temperature and precipitation). Alpine rock glaciers in the Pyrenees have been described by various authors, but to study them regionally has been a challenge since most of these studies are based on groundbased techniques. Two LiDAR campaigns (the first one completed in 2015 and the second one currently ongoing) performed by the Spanish Instituto Geográfico Nacional (IGN), up to now mainly used for land management purposes, present the perfect opportunity to combine both datasets and asses the movement of this periglacial features. The present work aims to show how by analysing two sets of open source LiDAR data (acquired in 2011 and 2016) and calculating elevation difference, the mapping of these features can be enhanced. The results will show how some recently proposed potentially active rock glaciers are showing vertical displacement and some other are not. Additionally, the displacement results of the Besiberri NW rock glacier are compared versus the results obtained by previous studies (1993-2003) by Chueca and Julian (2005) showing an increase on the vertical displacement of the central area interpreted as possible destabilization due to ice core ablation. To conclude, the climate monitoring implications, the limitations of the technique and further work to improve these results together with future applications are also discussed.


Introduction
Rock glaciers are typically defined as periglacial geomorphological features formed by angular rock debris bound together by ice, which can be massive or interstitial. Although the definition is widespread, there are still some significant discussions about their classification and origin, and whether they relate to purely glacial or periglacial environments since they are something "in between". Some authors argue that only debris flows bounded together by interstitial ice are to be considered as such, which will result in pure periglacial features (Haeberli, 1985;Barsch, 1992;Haeberli et al., 2006), meanwhile other authors defend that rock glaciers have a pure glacial origin, usually ice cored (Potter, 1972;Potter et al., 1998;Humlum, 1996). Aside from these discussions, what is broadly accepted is that these land morphologies represent a very solid permafrost indicator (Imhof, 1996;Barsch, 1996).
As stated before, rock glaciers are typically used as permafrost indicators, being permafrost defined as a ground thermal state in which, independently of if it is rock, soil or organic matter, it is found at 0 °C (or lower temperatures) for, at least, two consecutive years (van Everdigen, 1998). One of these scenarios is the alpine environments, which, independently from being located at lower latitudes, have specific locations which meet these criteria mainly due to the elevation and often aided by the orientation reducing the sun exposure during the warmest months of the year. The Pyrenees (Spain and France) are one of these cases (Serrano et al., 2009;Lugon et al., 2004;Fernandes et al., 2018).
The monitoring of Rock glaciers can be challenging due to its (generally) slow movement, which requires high precision measurements, and, up to date, have been mostly based on ground techniques. Differential GPS (dGPS) is an established technique, widely used for measuring the movement of rock glaciers (Schneider & Schneider 2001, Krainer & Mostler 2006. Remote sensing techniques offer a very powerful tool for analysing landforms, especially if combined by ground truth land-based surveys. Light Detection and Ranging (LiDAR) is a remote sensing technique which can be used from land-based and airborne platforms. Land-based LiDAR systems have gain popularity due to its high precision to monitor individual rock glaciers with repeated measurements at different times (Bauer et al., 2003) but have the disadvantage of the lack of wide area coverage, which makes them not fully suitable for regional studies. A recent example of the integration of land-based methodologies applied in The Pyrenees can be seen in Alonso-Gonzalez et al., 2018, which is applied to monitoring the Monte Perdido glacier, but not for rock glaciers.
On the other hand, traditional airborne techniques have been used in the past to analyse rock glacier movement using multitemporal data series, but mostly from optical image data comparison through photogrammetric techniques (Kääb et al. 1997, Kaufmann & Landstädter 2003, Kääb & Weber 2004, Roer et al., 2005. Airborne LiDAR data has been mostly used to assess morphological parameters of rock glaciers using Digital Elevation Models (DEM) obtained from LiDAR data (Janke, 2013). The added complexity in terms of cost and logistics (when compared to land-based LiDAR) has resulted in an frequent impediment for obtaining multitemporal datasets from the same platform to analyse rock glacier movement and inventorying. One of the few available papers which have achieved this is Bollmann et al., 2012, which focuses on the analysis of two rock glaciers.
Because of this, the present work aims to provide the first case study in the Spanish Pyrenees of quantifying vertical displacement. Also, how these results enhance the inventorying rock glaciers using airborne LiDAR multitemporal (two datasets from 5 years apart) for the Besiberri-Montardo and Punta Alta-Colomers area. For this, public data from the Spanish Instituto Geográfico Nacional (IGN) was used from two different LiDAR campaigns. Results are then compared against other authors, ground-based, vertical displacement measurements. Additionally, a further assessment of the potential rock glaciers identified by Ventura, 2016 within the study area is provided by analysing the vertical displacement and morphology of the movement. Furthermore, future work and regional applications in terms of inventorying and climatic impact analysis in rock glaciers for this workflow throughout all the Pyrenees are discussed, together with the limitations of the technique.

Regional Context
Location Figure 1: Location map showing the location of the currently inventoried rock glaciers included in this study. Capital letters represent measured and reported rock glaciers (green) and lowercase letters show potential rock glaciers (orange).
The study area is located in the eastern Spanish Pyrenees in the Lleida province. The bedrock on which the studied Rock glaciers lay is mostly made of granodiorites and monzogranites, with medium to coarse grains (Martin Parra et al., 2016). These Rock glaciers are located in the northern flanks of two mountain trends (NE-SW), also aligned with the Noguera Ribagorzana and Noguera de Tor river valleys, in which Caldas de Bohí is located. As a reference, the highest peaks found in the western lineation are Besiberri del Mig (3002 m), Besiberri N (3014 m) and Besiberri S (3030 m). On the eastern one, Pic de Contraix (2958 m) and Pic de la Pala Alta de Sarradé (2893m). The study area can be seen in Figure 1 together with the Rock glaciers included in this study (Table I).  Serrano et al., 2011, Chueca andJulian, 2005 and Ventura, 2016.

Rock glaciers
As mentioned in the introduction, extensive work has been done cataloguing rock glaciers in The Pyrenees. Serrano et al., 2011 synthesizes the work done and the rock glaciers identified. Chueca and Julian (2005) studied the Besiberri NW rock glaciers, which is included in this study. This tongue-shaped rock glacier is located in a glacial cirque, facing NW and is considered to be of glaciogenic origin (Chueca and Julian, 2005) related to the former Besiberri glacier plus additional talus material added by avalanches. García et al. (1998), by applying geophysical (geoelectric) surveys in the middle portion of the rock glacier, suggested the presence of a significant ice core (8-18 m) which is overlaid by a coarse-debris layer of approximately 1 to 1.5 m thick. Additionally, Chueca and Julian (2005) show how this ice is exposed at specific locations around the rock glacier.
By analysing geodetic data from a 10 year period, Chueca and Julian (2005) concluded that the horizontal flow rates vary between 3 to 24.6 cm yr -1 (8.72 to 13.35 cm yr -1 mean range) and subsidence estimations ranging from 1.5 to 13.8 cm yr -1 (mean ranges between 5 to 7.10 cm yr -1 ). In terms of horizontal velocity, an increase was noticed from the head to the toe area of the rock glacier, which they relate to an extensional flow dominating the rock glacier movement. Because of this, they also conclude that there is a likely deterioration of the ice-core.
Apart from these two Besiberri rock glaciers described with in situ observations and measurements, Ventura (2016) suggested, based on high-resolution orthoimages interpretation from 2D and 3D cartographic viewers a set of 69 potentially active rock glaciers not yet inventoried. Out of these 69, nine of them are within the study area and were assessed in this work. Ventura (2016) also classified these rock glaciers based on genetic and morphologic criteria; the genetic classification of these rock glaciers is also summarized in Table I.

Weather conditions
In terms of climatic conditions, the mean annual air temperature (MAAT) reported by Chueca and Julian, 2005, is ranging between -0.5 o C (at 3000m) to 6.8 o C at 1500 m and annual precipitation varying from 2600 mm (3000 m) to 1400 mm (1500 m). Additionally, and due to the lack of updated meteorological data in the area, the website Clima y Nieve Pirineos offers data from several different stations in and around the study area for the present year (up to May 2019) and the previous year (2018). To obtain a rough idea of temperature conditions in the study area, 2018 mean values from 11 stations at different elevations were used (Table II), and the temperature trend obtained was extrapolated to find an approximate elevation of the 0 o C isotherm, placing it at 2923.6m. It is important to consider that the real temperature at each location is most likely conditioned by local factors which cannot be taken into account with the available data, therefore, this estimated elevation for the 0 o C isotherm value, will be used as a qualitative reference to compare the Rock glaciers studied with. Furthermore, this obtained trend sets a temperature of -0.42 o C at 3000m in 2018, which is very close to what was described by Chueca and Julian, 2005. The correlation between the 2018 temperature data and the Rock glaciers elevation for the reported active in the Spanish Pyrenees plus some of the potential ones (Ventura, 2016) in or close to the study area in the Pyrenees can be seen in Figure 2.   Table II with a linear regression for Mean Annual Air Temperature showing the estimated 0 o C isotherm at 2923.6m.

LiDAR Data
The dataset used for this study is Open Access and has been acquired by the Spanish Instituto Geográfico Nacional (IGN), mainly for national land management purposes. The Light Detection and Ranging (LiDAR) program is composed of two different campaigns. The first one, which acquired data from all of the Spanish territory between 2008 and 2015, and a second one that started in 2016 and is still going on. The newest LiDAR datasets used in this study was acquired during one flight the 13 th of August 2016, meanwhile, the 2011 set was acquired in two separate flights on the 8 th of September and 2 nd of October 2011. All the technical specifications of these flights can be downloaded from the Centro Nacional de Información Geográfica (CNIG) website.
Both of these campaigns used an Airborne Laser Scanner (ALS) Leica ALS50 working with a minimum scanning frequency of 40Hz at a Field of View (FOV) of 50 o , minimum pulse frequency of 45kHz assuming a FOV of 50 o . The LiDAR system was used together with a Global Positioning System (GPS) and Inertial Measurement Unit (IMU). According to the website, the error specified for the 2011 campaign in the study area is of 0.3m (RMSE xy) and 0.4m (RMSE z); for the 2016 campaign is of 0.3m (RMSE xy) and 0.2m (RMSE z). In both cases, the point density is oscillating between 0.5-0.64 p/m 2 . In Figure 3 a summary of the coverage for both of these campaigns and the acquisition year is shown. The data was downloaded from the Download Center at the CNIG website and is supplied in compressed LAZ format. For each data sheet (a total of 34 sheets of 2x2km) two datasets are supplied which, according to the website, have gone through several Quality Control processes: (1) one with LiDAR Infrared (IR) elevation data, and (2) one in which points have been coloured in Red Green Blue (RGB) using simultaneous orthophotos from the same LiDAR flight.

Data processing and analysis
The first step of the process was to decompress the LAZ files into LAS format and, using ArcMap 10.3 (ESRI), convert these cloud points into Raster layers (using ETRS89, 31N UTM spatial reference). The RGB LiDAR data set was converted into a RGB raster layer and the IR elevation dataset was converted to both elevation and False RGB layers to end up with 3 images per year to work with. For this conversion, a binning approach was followed in which, for each output cell (of 1m cell size), was populated with the average of the points within each cell and filling the voids through triangulation using linear interpolation to determine the cell value. The result of this process is shown in Figure 4. Once the images are loaded and processed, an assessment of the general quality of the data was performed, which was considered of good quality. The 2011 set has very few data gaps, which are related to water bodies; meanwhile, the 2016 imagery shows more gaps, mostly related to water masses but also, in the north-western quadrant of the study area, to topographic highs. Additionally, some rectangular shapes with no data following a N-S orientation (these can be seen clearer in the IR False RGB image in Figure 4). Unfortunately, without the raw data from the LiDAR acquisition, it is not possible to understand what the reason for these gaps might be. Fortunately, most of the studied rock glaciers are not located in areas with data gaps except the most northern part of the tongue shape of Punta Harlé-2 and the NW lobe out of the two described by Serrano et al., 2011 in Besiberri N.
Secondarily, and since no raw and multiple return LiDAR data was available to understand the interaction of the pulsed laser light with snow-covered surfaces, a visual comparison between the RGB and IR False Colour was done. This comparison shows that textures and even features can be distinguished below the snow, which does not seem to have any relationship with snow cover ( Figure 5). This suggests that the images were acquired at a time where the snow cover was very local and probably very thin allowing the ALS obtain bare earth elevations correctly. Once a general quality assessment of both datasets has been achieved, the next step is subtracting the older (2011) elevation IR from the recent one (2016), resulting in an elevation difference product. Because of the data gaps described above and to enhance visually the small differences in elevation, the results were clipped and only values between -3.5m and 3.5m were used, setting the rest to Null Data. This threshold was selected based on the previous works (mainly by Serrano et al., 2006 and2010;Chueca andJulian, 2005 andGonzalez et al., 2011) in which the maximum vertical displacement per year reported is 22.5 cm yr -1 (Argualas Rock glacier, which is located 90km east from the study area) and 13.8 cm yr -1 for the Besiberri rock glacier within the study area. The threshold used is therefore considered suitable for imaging over 3 times the expected movement for five years. The resultant image can be seen in Figure 6 where one clear flow shape is highlighted (Besiberri NW rock glacier, in Figure 6a) and an example of noise related to vegetation and steeper slopes are shown also (Figure 6b). Ideally, the combination of Ground Control Points (GPC) and dGPS would provide an accurate estimation of the accuracy of the elevation difference calculated from the LiDAR datasets, but this information was not available. Because of this, although it is not ideal, the uncertainty assessment was calculated from the LiDAR cloud points individual error provided by the IGN, combined using a Root Mean Square Error (RMSE) estimation (1). Based on the RMS provided by the IGN on cloud point elevation, the RMSEzdiff diff estimated is 0.45m in absolute displacement. Considering this and that the background noise resulting from subtracting both raster layers is oscillating between -0.2m and -0.1m (Figure 7), the analysis of rock glaciers within this study will focus on analysing the elevation changes bigger than 0.5m. This background noise fits perfectly within the accuracy threshold specified in the technical specifications of the LiDAR acquisition. The fact that is not centred at cero could be related to slight calibration differences of the sensor for different flights. Additionally, to reduce the risk of interpreting incorrectly active rock glaciers, the morphology of the resultant deformation will also be taken into account together with the consistency of the deformation along the rock glacier shape identified from the RGB image.

Potential Active Rock Glaciers Identification
Based on the results obtained of elevation change (mostly subsidence), the two already identified active Rock glaciers (Southeastern lobe of Besiberri N and Besiberri NW) were reviewed showing, as expected, changes in elevation across their area ( Figure 8A and Figure 9E).  Regarding the potential active Rock glaciers defined by Ventura (2016), nine of them were reviewed, out of which four of them showed clear vertical displacement ( Figure 8B2, 8C, 8E and Figure 9D) with a consistent lobe or tongue morphology. Two did only show vertical displacement in very localized areas ( Figure 8B1 and Figure 9A and 9C), not generalized through the whole rock glacier seen in the RGB image. And two did not show significant displacement ( Figure 8D and Figure 9B), which could be related to no displacement or very little vertical movement hidden by the background noise. Additionally, this methodology allowed the addition of one new rock glacier to the inventory. In Figure 10, it is shown the displacement seen in one not included in the literature, except for Serrano et al., 2011, in which, in

Vertical Displacement
To provide an assessment of the magnitude of the displacement, the seven rock glaciers in which vertical displacement was identified were described more in detail. The vertical displacement along the profile shown in Figures 8 and 9 was averaged, but only considering the length within the profile that showed vertical displacement. This means that the length takes into account the general trend of vertical displacement (shown by the dashed line below the vertical displacement graph in Figures 8 and 9 for the rock glaciers with clear displacement identified). These parameters are summarized in Table III.

Integration of the results
Since most of the rock glaciers reviewed in this study have not been measured in terms of displacement, the only one to which the results obtained can be compared is Besiberri NW (Chueca and Julian, 2005). This comparison has been done by considering average displacement values for 3 transverse lines ( Figure 11) in two different ways: (1) averaging all the values within the transverse and (2) considering only the values at a regular spacing across the section, attempting to simulate how displacement measurements would have been done in the field in a similar way to Chueca and Julian (2005). These distances have been 20m for transverses A and B and 10m for transverse C, as shown in Table IV.  The results of this comparison (Table V) show a very good correlation with Chueca and Julian's (2005) results for transverse lines B and C. The section A-A' shows higher displacement (approximately double), which is interpreted as a destabilization of the central part of the rock glacier.  Chueca and Julian (2005), all the data along the transverse line and the measurements picked every 20m (A and B) and 10m in this study.
The results for sections B and C and its comparison with Chueca and Julian's (2005) suggest that this sections of the rock glacier have a velocity which has not significantly varied between the period of 1993 and 2003 measured by them and the 2011-2016 period studied here. This movement can be interpreted as regular creeping of the rock glacier downhill. On the other hand, the higher subsidence rate in transverse line A, suggest that there is an additional factor influencing the movement of this rock glacier activity. One of the possible explanations for this can be a quicker destabilization of the ice core described by Garcia et al. (1998) along the area surrounding section A, which can be also interpreted from Figure 11 by looking to the overall displacement morphology.
Additionally to the integration of the results obtained for the Besiberri NW rock glacier in previous literature, the resultant map can also be integrated with 3D visualizers such as Google Earth, resulting in very helpful graphical representation as shown in Figure 12.

Other possible periglacial features identified
Apart from the rock glacier diagnostic application discussed in this paper, this workflow can aid identify other features within this periglacial context. An example can be found in Figure 13, where two very clear scours which reach up to -2m elevation difference within the studied period. Also, some kind of levees are noticed with a maximum elevation increase of 0.8m. An additional increase of elevation at the bottom of these morphologies up to 1.2m maximum seems to be corresponding to sedimentation. These are currently interpreted as debris flows. This example shows that the workflow has huge potential, especially considering that the data is already available and open access or will be soon for the areas where the second LiDAR campaign is not completed yet. (4) Zoom on the 3D Google Earth TM Snapshot.

Implications in climatic monitoring
As widely known, rock glaciers act as a key diagnostic geomorphological feature when studying permafrost. Because of this, the monitoring of this features is highly important as it will provide a significant insight of how this climatic sensitive morphologies respond to environmental changes (temperature, precipitation and local conditions). They have the advantage, when compared to punctual and individual measurements, that these features respond to the overall interaction of different climatic indicators. Therefore, by monitoring changes in the behaviour of rock glaciers the health of permanent ice layer (independently from if it is massive or interstitial) within the rock glacier can be tracked.
The rough temperature assessment done in at the beginning of this paper shows that, when compared to what Chueca and Julian (2005) describe, there was an increase of 0.08 o C, which could be related to the different data sources (since theirs is not specified) and extrapolation techniques, but it could easily be related to the overall temperature increase in this alpine environment following worldwide heating trends widely known. This temperature increase could directly relate to the increase in subsidence rate described previously. This, could then be interpreted as a direct consequence of the ice core ablation due to generalized warming, which is something predicted by Chueca and Julian (2005) in their conclusions. Although the fact that this thinning of the rock glacier is not generalized down towards the toe is intriguing and worthy of further investigation since it might be related to diferential thickness of massive ice within the rock glacier.
Additionally, these type of Pyrenaic rock glaciers, located close to the 0 o C isotherm, are suitable for being very interesting natural 'sensors' for climatic change. Especially considering that the increase of velocity of rock glaciers in recent years have also been documented in the Alps (Bodin et al., 2018) and Northern Norway (Eriksen et al., 2018). Also, the fact that the multitemporal dataset used in this study is open access, easily available and will soon be covering all the Spanish mountain ranges, it can be widely used for assessing all the rock glaciers and understanding how they are responding to climate change in these past 5 years. If the LiDAR campaign is something recurrent that will be repeated every few years in the future, it will provide an excellent monitoring technique of climate change by considering not only the rock glaciers in this study, but all the available in the Spanish alpine environments (Pyrenees, Picos de Europa and Sierra Nevada) (Serrano et al., 2018).

Limitations of the technique
Probably one of the main disadvantages of the technique is that it does not allow the understanding of how the displacement rate varies from year to year, in other words, its temporal resolution. Since the assessment is based only on two elevation snapshots with five years' difference and the yearly displacement is averaged, it is impossible to know with the available data how this movement correlates with other climatic factors, like yearly temperature or precipitation changes for example.
Also, although the snow cover was assessed as one of the initial parts of the study, it is also a source of uncertainty since, as seen in Figure 14, some minor patches showing positive elevation changes corresponding to areas were snow or small puddles in small depressions are left over by meltwater. This response is very similar to the one given in water bodies. Fortunately for this study, this only occurs in very specific locations, and as seen in Figure 14, it does not compromise the results obtained since the subsidence patterns do not seem to have any correlation with these snow patches. Furthermore, the error estimated for this workflow, when compared to the yearly displacement calculated, makes this workflow considering the airborne LiDAR datasets not fully suitable for short term elevation displacement comparison. This makes the methodology a very interesting tool for regional analysis of these type of landforms, but for more detailed, yearly studies, complementary tools and methods should be explored.
Apart from all of this limitations, the fact that the data does not imply an extra cost and that the methodology is fairly simple and straight forward, the advantages compensate significantly the disadvantages for extrapolating this workflow into other areas and maybe applications, like geomorphologic cartography.

Further work to be done
Since this methodology has been proven effective to aid in the inventorying of rock glaciers and additional geomorphological features, it could be easily extrapolated to the rest of the Spanish mountain ranges. This implies that, when the full second LiDAR campaign is available, it can aid in identifying these features and therefore, help to update the alpine permafrost cartography proposed by Serrano et al., 2009. Also, as mentioned before, to keep track of the climatic impact of these features, since they are very sensitive to temperature and precipitation conditions.
Additionally, the vertical displacements here described are suitable for analysis using Interferometry from Synthetic Aperture Radar (InSAR) methodologies from Sentinel-1 or TerraSAR-X data, which will aid in the identification of moving features, increasing the confidence and accuracy in the mapping. Additionally, since these satellites cover the same area weekly and biweekly, it can be very interesting to analyse how the movement of these features varies from year to year and even throughout one same summer-autumn season (since snow cover will affect the measurements) and compare these results against temperature and precipitation yearly variations.
Another airborne technique very interesting is the use of commercial drones or Unmanned Aerial Vehicles (UAV) and photogrammetric techniques, which can aid in the detailed study of these features individually. With these methods, in combination with dGPS Ground Control Points (GCP) and Real-Time Kinetics (RTK) or Post-Processing Kinetics (PPK), the accuracy can be down to a few centimetres. With repeated surveys every year, the evolution of these rock glaciers can be very nicely mapped and studied, since the 3D models obtained can be analysed through multiple spatial analysis tools available. This will also provide an assessment of horizontal displacement and the real velocity of the feature. This will result in a very powerful perspective, complementary to the LiDAR data regional analysis.

Conclusions
Through the analysis of LiDAR data, publicly available, intended mainly for land organization purposes, this study has proven how it can aid in the identification of permafrost-related geomorphological features. This adds significant value to the dataset which, up to now, it is believed it has not been used for this purpose.
The analysis has helped to update the rock glacier inventory in the Besiberri-Montardo and Punta Alta-Colomers area providing an assessment on the movement of potential rock glaciers proposed by Ventura (2016). From the nine potential rock glaciers proposed, this study proved that four of them show generalized clear vertical displacement with tongue or lobe morphologies, and therefore, are correctly categorized as active rock glaciers. Two show minimum movement in very localized areas of the rock glacier and other two does not show any movement above the background noise described in the methodology, which can mean that no movement is taking place or that it is below the resolution of the method.
Additionally, the Besiberri NW, studied in situ by Chueca and Julian (2005), showed very similar average displacement across two out of the three sections compared. The one differing from their results, with an increased subsidence rate at the central part of the rock glacier, is suggested to be related to the ablation of the ice core described in previous studies.
To conclude, the possibilities of this type of workflow are highlighted together with alternative methodologies which can be complementary and can aid in the identification and monitoring of these permafrost-related features. The application of this monitoring concerning the current climatic variation was also discussed, especially considering that the data used in this study will be available for all the Spanish territory soon.