Combination of Aerial, Satellite, and UAV Photogrammetry for Quantifying Rock Glacier Kinematics in the Dry Andes of Chile (30°S) Since the 1950s

The diachronic analysis of aerial and satellite imagery, uncrewed aerial vehicle (UAV) and in situ surveys obtained between 1956 and 2019 are employed to analyse landform surface kinematics for the Tapado site located in the Dry Andes of Chile. A feature tracking procedure was used between series of orthorectified and co-registered images to calculate surface velocities on several ice-debris landforms, including rock glaciers and debris-covered glaciers. For the active rock glaciers, the results exhibit typical viscous flow, though local destabilisation process seems to occur, increased velocities since 2000 (>1 m/yr) and terminus advance. Nevertheless, the debris-covered glaciers displays heterogeneous spatial patterns of surface velocities, together with collapse (downwasting) associated with the development of thermokarst depressions and supraglacial ponds. Our findings show that surface kinematics and multitemporal observations derived from different sensors are valuable tools for differentiating between glacial and periglacial features. The pluri-decadal time series since 1956 constitute a unique dataset for documenting the surface kinematics of creeping mountain permafrost in the Southern Hemisphere. The approach developed in this work offers a way forward to reconstruct the recent behaviour of glacial and periglacial features in the Andes, where archival aerial photographs are available but have not previously been processed rigorously to obtain an accurate assessment of landform kinematics.

The diachronic analysis of aerial and satellite imagery, uncrewed aerial vehicle (UAV) and in situ surveys obtained between 1956 and 2019 are employed to analyse landform surface kinematics for the Tapado site located in the Dry Andes of Chile. A feature tracking procedure was used between series of orthorectified and co-registered images to calculate surface velocities on several ice-debris landforms, including rock glaciers and debris-covered glaciers. For the active rock glaciers, the results exhibit typical viscous flow, though local destabilisation process seems to occur, increased velocities since 2000 (>1 m/yr) and terminus advance. Nevertheless, the debris-covered glaciers displays heterogeneous spatial patterns of surface velocities, together with collapse (downwasting) associated with the development of thermokarst depressions and supraglacial ponds. Our findings show that surface kinematics and multitemporal observations derived from different sensors are valuable tools for differentiating between glacial and periglacial features. The pluri-decadal time series since 1956 constitute a unique dataset for documenting the surface kinematics of creeping mountain permafrost in the Southern Hemisphere. The approach developed in this work offers a way forward to reconstruct the recent behaviour of glacial and periglacial features in the Andes, where archival aerial photographs are available but have not previously been processed rigorously to obtain an accurate assessment of landform kinematics.

INTRODUCTION
Cryosphere components in high-mountain regions such as glaciers, snow cover and permafrost are sensitive to environmental changes such as rising temperatures because of their proximity to melting or thawing conditions (Haeberli and Beniston, 1998;Barry and Gan, 2011;Hock et al., 2019). Permafrost is invisible for the layperson because it is a thermal phenomenon (Dobiński, 2011); however, direct temperature measurements (Noetzli et al., 2021), geophysical measurements (Mollaret et al., 2019), geomorphological evidence (Barsch, 1996) and modelling approaches (Gruber, 2012) can reveal permafrost conditions in diverse mountain or polar environments. In the Dry Andes (17°30´-35°S, sensu Lliboutry, 1998), rock glaciers are regarded as the unequivocal expression of the creep of mountain permafrost, and it is assumed that they play an essential role as water reservoirs due to their prevalent occurrence in the high Andean landscape (Corte, 1976;Schrott, 1996;Azócar and Brenning, 2010;Perucca and Esper Angillieri, 2011;Schaffer et al., 2019). In this mountain region, rock glaciers, push moraines and debris-covered glaciers assemblages represent cases of intricate permafrost-glacier interactions (Trombotto et al., 1997;Bodin et al., 2010). In addition, recent studies performed in the Dry Andes have highlighted the remarkable development of several transitional landforms from glacial to periglacial domains Kinnard, 2015, 2017). These studies stress glacial and periglacial processes' importance and encourage more in-depth investigations to characterise this hybrid cryospheric landscape. Nevertheless, and due to this complexity, glaciers and periglacial landforms are not distinguished from one another in current country-wide glacier inventories at the Chilean and Argentinian share of the Andes (Barcaza et al., 2017;Zalazar et al., 2020). Environmental laws designed to protect surface and subsurface perennial ice reservoirs (sensu stricto glaciers and rock glaciers) have been implemented in Argentina, and are currently in discussion in Chile, which means that their correct characterisation and monitoring have become a highly relevant environmental matter (Brenning, 2008;Brenning and Azócar, 2010;Kronenberg, 2013;Wainstein et al., 2020).
In the Andes, one of the earliest rock glacier monitoring activities was performed in 1970 at the Pedregoso rock glacier in the Valparaiso Region of Chile (32°S). Using a triangulation land survey technique, Marangunic (1976) measured 23 points during the late summer of 1970, estimating velocities between 0.2 and 1.28 cm per day at the rock glacier root and the central part, respectively. Later on, Valenzuela (2004), reported the kinematics of South Infiernillo Alto rock glacier near the Los Bronces mining site (33°S). In this particular case, the rock glacier reached up to 25 m/yr in 1995 after the deposit of 14 m tons of mining waste on the rock glacier's central part. Terrestrial geodetic surveys have been sporadically implemented at the Tapado rock glaciers since 2009 (Dirección General de Aguas [DGA], 2010), between 2012 and 2016 at the Varas rock glacier (23°S) in the Norwest Argentinean Andes (Martini et al., 2017) and between 2013 and 2016 at the Morenas Coloradas rock glacier (32°S) in the Central Argentinean Andes (Trombotto-Liaudat and Bottegal, 2020). More recently, Blöthe et al. (2020) provided the surface velocities of 244 rock glaciers between 2010 and 2018 in the Cordón del Plata Range (Central Argentinean Andes) using highresolution satellite images and cross-correlation matching techniques. Using a combination of UAV and geophysical surveys, Halla et al. (2021) investigated the surface kinematics, ice content and interannual water storage at the Dos Lenguas rock glacier, approximately 15 km to the southeast of the Tapado area. Despite the recent upsurge of rock glacier studies in the Andes, the lack of long-term systematic measurements of these landforms has precluded a better understanding of their response to changes in environmental conditions (Masiokas et al., 2020) and more appropriate classifications of rock glaciers according to their degree of activity (i.e., the traditional distinction between active and inactive landforms, sensu Barsch, 1996).
Several studies have highlighted the coherent nature of rock glacier surface deformation over time, indicating a large-scale stress transmission mechanism on the ice-debris mixtures under permafrost conditions (Wahrhaftig and Cox, 1959;Kääb, 2002;Haeberli et al., 2006). This coherent surface deformation is mainly the result of permafrost creep processes occurring inside the rock glacier body, like the plastic deformation of frozen sediments supersaturated with ice and the localised deformation at the shear horizon (Wagner, 1992;Haeberli, 2000;Arenson et al., 2002;Cicoira et al., 2021). As a result of these intrinsic rock glacier dynamics, changes in surface morphology and kinematics fields are dissimilar between rock glaciers and other ice-debris complexes like debris-covered glaciers (Berthling, 2011;Bosson and Lambiel, 2016). However, differentiating among debris-covered glaciers, glacial and periglacial landforms from solely monotemporal remote sensing observations is non-trivial (see, for instance, the classification scheme proposed by Janke et al., 2015), and the interpretation of composite geomorphological processes is potentially ambiguous. Therefore, multitemporal observations based on archival or contemporary high-resolution imagery are more suitable to describe complex geomorphological processes in high mountain environments (Roer and Nyenhuis, 2007;Micheletti et al., 2015;Groh and Blöthe, 2019;Hu et al., 2021;Kääb et al., 2021).
Conventionally, rock glacier inventories and classification schemes are based on the presumed degree of activity by the visual inspection of optical imagery alone (Brardinoni et al., 2019). Furthermore, deliberations on rock glacier activity based only on surface morphology indicators (see Imhof, 1996;Burger et al., 1999) should be taken with precaution, as controversy and confusion have intensified from morphological centred approaches on complex landform mapping (Berthling, 2011). Currently, kinematical data obtained from optical and radar remote sensing datasets (Delaloye et al., 2007;Barboux et al., 2014;Villarroel et al., 2018;Hu et al., 2021;Kääb et al., 2021) has been recognised to be instrumental in distinguishing among different types of slope movements in the high mountain periglacial environments. Typically, rock glacier kinematics can be evaluated from remotely sensed datasets, either airborne or space-borne, obtained for at least two different points in time to calculate significant displacements above a certain threshold (Kääb, 2005).
In the present study, we aim to analyse spatial-temporal changes of surface kinematics over rock glaciers and ice-debris complexes in the Dry Andes by exploiting an assemblage of remote sensing datasets spanning 63 years . This work supports a geomorphological mapping proposal based on multitemporal remote sensing techniques and geomorphic interpretations informed by observations on the evolution of landform kinematics. Furthermore, we provide a systematic description of the development of rock glacier kinematics over distinct types of landforms since 1956. Our work is embedded into the efforts of the recently founded International Permafrost Association (IPA) Action Group (RGIK, 2021) on rock glacier inventories and kinematics by deciphering long-term surface velocities.

STUDY SITE AND PREVIOUS WORK
The study site is located in the La Laguna basin at the headwaters of the Elqui River catchment in the Coquimbo Region of Chile (30°S, 69°W, see Figure 1). The regional climate is characterised by semiarid conditions and is influenced mainly by the extratropical precipitation systems (Garreaud, 2007), although nonnegligible convective activity is occasionally present during the austral summer months (Vuille and Keimig, 2004). The annual precipitation amounts and corresponding snowpack processes exhibit substantial variability associated with El Niño-Southern Oscillation (ENSO, e.g., Montecinos and Aceituno, 2003;Réveillet et al., 2020), with above (below) -average precipitation during El Niño (La Niña) events. At the nearest long-term meteorological station (La Laguna, 3,160 m a.s.l.), the mean annual air temperature (MAAT) was 8.15°C in 1976-2019; and the mean annual precipitation was 156 mm in 1964-2019. The MAAT at the La Laguna station experienced a warming trend of 0.17°C per decade . Furthermore, a sequence of dry years has been observed since 2010 in this region, with a precipitation deficit between 20 and 40% (Garreaud et al., 2020). Short-term meteorological data from an automatic weather station (AWS) at 4,306 m a.s.l. (Figure 1) indicates a MAAT of -0.6°C between 2014 and 2020 (CEAZA, 2021). According to a recent permafrost distribution model developed by Azócar et al. (2017), permafrost is scattered between 3,900 and 4,500 m a.s.l. and more prevalent above 4,500 m a.s.l.
The Tapado area ( Figure 1) includes the Cerro del Tapado (5,536 m a.s.l.), which is covered with the largest glacier in the region (Barcaza et al., 2017). The Tapado Glacier has drawn much research attention in recent years due to its relatively easy accessibility, large area and unique position south of the socalled "South American Arid Diagonal" (Nicholson et al., 2016;Pourrier et al., 2014;Sinclair and MacDonell, 2016). The analysis of a 36 m long ice core and borehole temperatures near the summit retrieved in February 1999 revealed a large interannual variability of net accumulation during the 20th century and temperatures down to -11°C at the ice-bedrock interface (Ginot et al., 2006). In connection with the Tapado Glacier, a complex glacio-geomorphologic assemblage (here termed Tapado complex) of rock glaciers, moraines, and a debris-covered glacier is present . Next to the Tapado complex, Las Tolas and El Cachito rock glaciers ( Figure 1) display clear signs of ongoing activity such as transversal ridges and furrows and steep lateral-frontal talus ( Figure 2). The dominant surface material in this area is mainly a mix of pebbles and cobbles derived from the Pastos Blancos and Doña Ana formations, whose lithologies are mainly silicic volcanic rocks (Strauch et al., 2006). According to their predominant clast size, the rock glaciers are classified as pebbly rock glaciers (sensu Ikeda and Matsuoka, 2006), which is also the dominant typology in this part of the Andes (Falaschi et al., 2014). Earlier reconnaissance work done by Paskoff (1967) detailed the first geomorphological map of the Tapado site as well as some initial insights into the complex interplay between glacial and periglacial landforms. Later on, Milana and Güell (2008) performed seismic refraction surveys on the lower Tapado complex and Las Tolas rock glacier, inferring active layer thicknesses of 4 and 11 m, respectively. From these findings, they postulated that the Tapado rock glacier was derived from the debris-covered glacier and labelled it of "glacigenic" origin based on the higher ice content detected by their geophysical survey. Monnier et al. (2014) performed more extensive geomorphological and geophysical measurements, concluding that the debris-covered glacier overlapped the Tapado rock glaciers after a noticeable glacier advance, which probably occurred during the Little Ice Age (LIA) climate episode (Espizua and Pitte, 2009). Recent studies have also focused on the hydrological mechanisms operating through the different cryospheric components (Pourrier et al., 2014) and their hydrological relevance (Schaffer et al., 2019).

MATERIAL AND METHODS
The methodological approach provides an integrated kinematic quantification and geomorphological analysis using various multi-temporal remote sensing datasets ( Table 1). Distributed kinematic changes were quantified at different time intervals within the period 1956-2019. Products of the kinematic quantification were subject to a quantitative evaluation and a detailed geomorphic interpretation. The ensemble of orthoimages was projected to the Universal Transverse Mercator (UTM) 19S zone and WGS-84 Datum to provide compatible datasets. Elevation values were recorded in heights above the WGS-84 ellipsoid.

Airborne Laser Scanning Dataset
The airborne laser scanning (ALS) survey was performed by Digimapas company on 18 April 2015, using a Trimble Harrier 68i laser scanning system onboard a helicopter, flying at 640 m above the ground. The airborne system was embedded with a metric camera, an inertial measurement unit (IMU) and GNSS devices. The postprocessing (trajectories, point cloud filtering, The lack of surface contrast caused by a thin layer of seasonal snow excluded the use of the accompanying orthomosaic for the ALS 2015 dataset for further kinematic measurements. However, this dataset proved to be extremely useful during the extraction of photogrammetric ground control on assumed stable areas outside the Tapado complex. Following the approach presented by James et al. (2006), several ground control points (GCPs) were extracted from the ALS dataset based on two main principles: 1) correctly and easily recognisable on several images; and 2) located on stable zones during each period. These requirements were typically met by medium-sized boulders on flat areas and rock outcrops near mountain ridges.

Geodetic Measurements
The surface kinematics of the Tapado rock glacier have been surveyed on a network of 46 points since 2009 using differential GNSS (dGNSS) system (Dirección General de Aguas [DGA], 2010). A temporary base station on a flat exposure of bedrock, near the rock glacier, was deployed during each surveying campaign with a logging interval of 1-2 s. Additionally, around three control points located on stable terrain outside the rock glaciers were measured during each campaign as a quality check and to ensure the stability of the base station. During the 2018 and 2019 field surveys, this control network was augmented to include between 14 and 17 GCPs to be used during the processing of the UAV surveys. The base station coordinates were fixed using the Trimble Center Point RTX postprocessing service, and the differential postprocessing of the GNSS raw data between this base and the rover GNSS antenna was conducted using Trimble Business Center (TBC, V.4) surveying software. The stability of the base station and the coordinate system is discussed in section 5.3. The reported average horizontal and vertical precisions (95%) were 0.02 and 0.04 m, respectively.

Archival Aerial Photographs
Archival aerial photographs from the 1956 HYCON aerial survey were provided by the Instituto Geográfico Militar (IGM) of Chile, and the 1978 CHILE60 and 2000 GEOTEC aerial surveys were provided by the Servicio Aerofotogramétrico (SAF) of Chile. In the high Chilean Andes, and until recently, the general approach for analysing raw archival aerial photographs sourced from either the IGM or SAF was their treatment by georeferencing using the available regular 1:50000 scale cartography or recent datasets (Bown et al., 2008;Rabatel et al., 2011;Ruiz Pereira and Veettil, 2019; among others). However, georeferencing itself cannot compensate for the distortions induced by the camera optics and the extreme relief encountered in such mountain regions (Mikhail et al., 2001). Therefore, more robust processing based on reconstructing the interior and exterior image orientations of historical aerial photographs can be sought using modern digital photogrammetric techniques and high-quality ground control (see Farías-Barahona et al., 2019. With the reconstruction of three-dimensional terrain data through stereo images, input imagery can be transformed into geometrically corrected images (i.e. orthoimages and orthomosaics). We applied this methodology using Geomatica Banff photogrammetric software.
Aerial photographs were scanned at 1,200 dpi and 8-bit/pixel using a photogrammetric scanner and provided with their corresponding calibration certificates, except the 1956 photographs for which the camera calibration certificate was not available. Therefore, we only retrieved the calibrate focal length from the 1956 photographs' data strip and estimated the rest of the interior orientation parameters such as principal points and radial distortion parameters through a self-calibration adjustment implemented in Geomatica Banff photogrammetric software. Such adjustment requires more GCPs than when calibration certificates are available. A bundle block adjustment ( Table 2) was performed based on several GCPs extracted from the 2015 ALS dataset (point cloud and orthomosaic) and automatically generated Tie Points (TPs) followed by a subsequent aerotriangulation, including adjusted camera parameters. Using a Semi-Global Matching algorithm, DEMs were extracted from stereo images at a resolution of 1 m × 1 m. Orthoimages were generated at the same resolution and extent (i.e. concurrent) using their corresponding DEMs.

High-Resolution Satellite Imagery
Two stereo panchromatic images from the GeoEye-1 and two tristereos Pléiades satellites were available for this study ( Table 1). GeoEye-1 and Pléiades images were provided in level 1B and The area covered by the UAV surveys is indicated in Figure 1.
Frontiers in Remote Sensing | www.frontiersin.org November 2021 | Volume 2 | Article 784015 coded over 11-bit/pixel and 12-bit/pixel for each sensor, respectively. The original orientation data obtained from the rational polynomial coefficient (RPC) was refined using 8 GCPs and between 65 and 85 automatically extracted TPs. The GCPs were sourced from the 2015 ALS dataset similarly to the aerial photographs. Table 2 shows the accuracy of the GCPs and TPs for each bundle block adjustment. Using the OrthoEngine module in Geomatica Banff, the processing of the satellites images follows the same preparation as the aerial photographs to generate high-resolution orthoimages and DEMs for each year.

Uncrewed Aerial Vehicle Surveys
The UAV surveys were planned to acquire recent high-resolution images of the northern frontal lobe of the Tapado complex (  4 Pro DJI devices (Table 1), and their planning and execution were achieved employing the PIX4Dcapture flight planning app. The front and side overlaps were set to 80 and 70%, respectively, using a single grid flight plan. The nominal altitude above the ground during both surveys was kept around 45 m by flying nearly parallel to the slope, and the final ground sampling distances (GSD) were below 2 cm. The total area covered by both UAV surveys reached 0.2 km 2 . The UAV images were processed using PIX4Dmapper pro version 4.4 software through a prescribed Structure-from-Motion (SfM) photogrammetry approach to derive 3D point clouds, DEMs and orthomosaics (Carrivick et al., 2016). Noise filtering and medium surface smoothing were activated in the software processing options. A bundle block orientation of each survey was achieved using 14 (2018) and 17 (2019) GCPs distributed on the rock glacier surface (Figure 3). The resulting orthomosaics and DEMs were sampled at a resolution of 2 cm × 2 cm. High-resolution elevation changes were computed by the DEMs of differences (DoD) method, using the Geomorphic Change Detection (GCD), version 7 software (Wheaton et al., 2010). Between 16 (2018) and 11 (2019) checkpoints (CPs) were strictly employed to independently gauge the quality of their corresponding DEM (Figure 3). The standard deviation of the elevation differences between the reference CPs and DEM was used to calculate the minimum level of detection (LoD), which determines the threshold between significant and non-significant elevation change (at 95% confidence interval).

Feature Tracking
Each orthoimage pair was analysed at their respective coarsest resolutions (see Table 1). Surface feature tracking was obtained by a normalised cross-correlation procedure using CIAS software (Kääb and Vollmer, 2000;Debella-Gilo and Kääb, 2011). The dimensions of the reference and search windows were first defined based on the orthoimages resolution, quality and time interval considered. Typically, smaller (larger) reference and search windows were employed for short (long) periods. The resulting sizes of the reference (search) windows are 64 × 64 (128 × 128) for aerial images (early periods), 25 × 25 (50 × 50) for satellite images (recent periods) and 50 × 50 (200 × 200) for UAV-derived orthomosaics. To derive comparable surface velocities, conjugated (corresponding) points in two orthoimages were spaced within the predefined boundaries based on an invariable 10 m sampling grid (i.e. eulerian framework). Older aerial photographs presented different image quality due to dust, dirt and scratches on the film during the scanning process, reducing the quality of the feature matching processing. On the other hand, modern high-resolution satellite imagery presented better conditions for feature tracking due to their higher radiometric resolution of 11-12 bit/pixel and better signal-to-noise ratio with respect to the archival aerial photographs scanned at 8-bit/pixel (Poli et al., 2015).
The accuracy of the kinematic time series was assessed by stable features such as boulders and outcrops that were considered representative of the study landforms ( Figure 4). This assessment is done by matching features on stable ground, computing an apparent x-y shift, scale and rotation values between the two images, and then the final residuals after the adjustments (Kääb, 2021). The x and y standard deviations (σ x and σ y ) of the residuals were used to calculate the uncertainty associate with the horizontal displacements (σ l ) following the rationale presented by Redpath et al. (2013) and evaluated using the LoD with a confidence limit of 90% (i.e. 1.64×σ l ). The resulting LoD values were converted to m/yr based on the time difference between subsequent datasets ( Table 1).
Different conditions were applied to filter and exclude poorly correlated points and those displaying anomalous displacement or directions. Depending on the intrinsic quality of orthoimages pairs (e.g. image 1 and image 2), reliable displacement values were retrieved with correlation coefficients higher than 0.4-0.6 (Wangensteen et al., 2006;Kääb et al., 2021). Zones with significant changes in the surface texture such as thermokarst depressions and ponds, fluvioglacial fans, snow patches and exposed glacier ice hampered the feature tracking procedure. Therefore, these zones were excluded from the feature tracking analyses. Time series of surface kinematic were derived using the mean value within manually drawn kinematic zones, defined by an almost coherent movement and the availability of high-quality correlation points during the analysed periods.

Ground and Air Temperature Data
Ground temperatures from two shallow boreholes were available at the northern frontal lobe of the Tapado complex (Figure 1). The setup consists of five temperature data loggers (E348-S-TMB-M006 HOBO) placed at depths between 0.2 and 2.5 m (Dirección General de Aguas [DGA], 2010). Recordings are stored in a centralised datalogger (HOBO U-30) with a logging interval of 30 min. Air temperature data was also measured on each borehole site, but the recordings were quite discontinuous. Long term mereological data was sourced from

Evaluation of the Displacement Accuracy
The horizontal accuracy assessment between orthoimage pairs was achieved by evaluating the apparent displacement of the stable points outside the studied areas ( Figure 3 and Table 3). As expected, the older orthoimage pair presents the larger LoD associated with the more significant errors during the bundle block adjustments ( Table 2) and the relatively low quality of the original scanned film. Additionally, the 1956 aerial survey is the only one where the internal calibration parameters were reconstructed using a self-calibration technique. Nevertheless, this large LoD is somehow compensated by the significant time span between the early aerial photography acquisitions (nearly 22 years). Smaller LoDs are associated with the orthoimages pairs from the recent satellite acquisitions (GeoEye-1 and Pléiades) and the UAV-derived orthomosaics. The standard deviations associated with CIAS analysis were approximately equal in the y and x directions for most orthoimages pairs, revealing a small anisotropy in the resulting displacement uncertainties.   textural changes. Surface velocities below the LoD were concentrated in the recently deglaciated glacier forefield and near the unit's margins. Velocities of more than 1.5 m/yr were detected in the upper debris-covered glacier; and the central and lower sections of Las Tolas and El Cachito rock glaciers, respectively. The velocity profiles of the three most active rock glacier units can be seen in Figure 6. The fastest rock glacier part was up to 1.8 m/yr near the front of the El Cachito rock glacier. The Las Tolas rock glacier presented the fastest section (1.6 m/yr) around the middle approximately halfway between the rooting zone and the terminus. On the Tapado 1 rock glacier (T1rg), the increase in velocity along the centre profile was neither uniform nor linear but indicated at least two rapid sections ( Figure 6).

Surface Kinematics
Time series of surface kinematics were compiled for two rock glacier units at the Tapado complex (T1rg and T2rg), Las Tolas (LTrg) and El Cachito (ECrg) rock glaciers (Figure 7). However, between 1956 and 1978, results from CIAS analysis at the El Cachito rock glacier showed values below the LoD; therefore, no kinematic values were reported. Figure 7 also includes the estimated MAAT at 4,500 m a.s.l. since 1976 and the recent in-situ kinematic monitoring by GNSS surveys at the Tapado 1 rock glacier (2010-2019). Excluding the Tapado 2 rock glacier, the group of landforms display a nearly synchronous kinematic behaviour. The Tapado 2 rock glacier displayed nearly constant velocities and low activity (<0.3 m/yr, see Figure 7). Surface kinematics indicate relatively low velocities in the 1950-2000 period, followed by a rapid acceleration since the 2000s. Over the period 2010-2012, the surface kinematics of the El Cachito, Tapado 1 and Tapado 2 rock glaciers experienced a slight deceleration, which lasted until 2015 for Tapado 1. The Las Tolas rock glacier was the only landform that did not decelerate during the entire study period.
The 63 years glacier retreat and periglacial processes can be observed from the dynamic visualisation of the sequential orthoimages at the video supplement (Supplementary video S1). On the one hand, this dynamic visualisation also shows the chaotic downwasting structure and the development of thermokarst depressions and ponds at the debris-covered glacier. On the other hand, the viscous flow appearance of

High Spatial Resolution Surface Changes
The spatial surface behaviour of the Tapado 1 rock glacier is analysed based on the recent UAV-derived DEMs and orthomosaics ( Figure 8). The delineation of the surface between the central and frontal parts and lateral talus (Figure 8) was distinguished by the distinctive elevation change patterns (see Halla et al., 2021 for a similar rock glacier division). From the high-resolution DoD, the advection of transverse ridges is seen through the alternating positive and negative elevation changes at the rock glacier central tongue. Here, surface rising (sinking) patterns are found in the front (rear) of individual ridges. The amplitude associated with the surface rising and sinking can reach between 0.5-1 m, whereas the wavelength of the ridge-and-furrow morphology is around 20 m. The widespread positive elevation changes at the frontal and lateral slopes are associated with the advance of the rock glacier and reworking processes. Negative elevation changes (up to -1.5 m) are associated with some rockfalls and gullying patterns just below the rock glacier front. The recent changes of the Tapado 1 rock glacier can be observed from the dynamic visualisation of the sequential highresolution hillshades at the video supplement (Supplementary video S2).
The surface velocities at the Tapado 1 rock glacier derived from recent satellite images (2014-2019) and UAV (2018-2019) display a similar structure (Figures 5 and 8). Surface velocities over 1.25 m/yr were measured near the rock glacier front. Near the lateral margins, there is a sharp decrease in the observed velocities. Surface velocities remained high near the rock glacier front. Due to the steep topography of the frontal and lateral talus (slopes angles between 32 and 39°) and the significant texture changes, no reliable velocities values were obtained on this zone. Nevertheless, visual inspection of the high-resolution orthomosaics indicates that the surface changes are associated with erosive processes such as falling boulders, small debris slides and localised gullying.

Landform Interpretation
The following landform interpretation is mainly based on the recent guidelines published by the IPA Action Group on rock glacier inventories and kinematics (RGIK, 2021). Some subjectivity is still inherent in the difficulty of conciliating kinematics zones with fuzzy geomorphological boundaries.

Tapado Complex
Following the recent definitions of spatial connection of the rock glacier to the upslope unit, the Tapado 1-5 rock glacier units are glacier-connected ( Figure 9) because they interact with the debriscovered glacier and push moraine complexes. Independently from the origin of the ground ice, glaciers can supply large amounts of meltwater and ice to proglacial permafrost landforms. As the case for Tapado 6 rock glacier, this unit is classified as glacier forefieldconnected, which probably was in contact with the Tapado glacier during the LIA advance. This recently exposed glacier forefield is the heritage of the former Tapado glacier dynamics, and its thermal regime (i.e. polythermal glacier) probably explains the presence of several fluted moraines, fluvioglacial fans and dry gullies. Due to the diffuse margin between this rock glacier unit and the glacier forefield, the Tapado 6 rock glacier unit was previously interpreted as part of an upper lateral moraine complex . However, the surface velocities up to 0.25 m/yr and the well-developed frontal and lateral talus indicate an active rock glacier unit, though the development of the landform remains difficult to determine when considering the evolution of the glacier.
Large thermokarst depressions might indicate the boundary between the Tapado debris-covered glacier and the push moraines ( Figure 9). These push moraines (glaciotectonized frozen sediments) are associated with the former glacier advance during the LIA. Following glacier retreat and downwasting, these push moraines have developed a particular back creeping process toward the former glacier body (Haeberli, 1979;Kneisel and Kääb,  2007), which is shown in Figure 9. Despite the short transport route of fine-grained sediments in Tapado glacier forefields, the size of this push moraine formation is quite exceptional and may indicate the incorporation of a preceding ice-debris landform during its development.

Las Tolas and El Cachito rock glaciers
The Las Tolas rock glacier is a talus-connected landform. The topographical sequence goes from a massive headwall of an icefree cirque with several coalescent talus slopes, a distinguished body of permafrost creep with the typical ridge-and-furrow morphology, and ending with a trident-like form ( Figure 5).
Considering the particular topographical setting, the El Cachito rock glacier was catalogued as a debris-mantled slope-connected rock glacier. Multitemporal observations and kinematic data suggest the debris are provided from the in-situ weathering process founding the upslope geomorphological unit and subsequently mobilising downslope to convene into a wellorganised expression of creeping mountain permafrost (Supplementary video S1).

1 Potential Sources of Uncertainty
The study area is situated on a tectonically active mountain range due to the Nazca plate's subduction beneath the South American one. Several earthquakes have hit the area during the last decades, generating large coseismic and postseismic displacements (Ruiz and Madariaga, 2018). Remarkably, the September 2015 Mw 8.3 Illapel earthquake generated coseismic displacements up to 2.13 m across the whole region (Barnhart et al., 2016). The study site is located 235 km NE from the earthquake epicentre. Postprocessing values of the Tapado GNSS base station before (2014) and after (2018) the Illapel earthquake indicated a westward motion of 0.2 m, which is coherent with the data from the nearest earthquake monitoring GPS station (Shrivastava et al., 2016). Together with the continuous mountain uplift, such displacements have generated various horizontal and vertical shifts in the global coordinate system among fixed points, which were assumed to remain stable. However, it should be remarked that this study involved using a standard set of GCP and CP (extracted from the ALS and GNSS surveys before September 2015) during the orthorectification and coregistration procedures. This procedure ensured that orthorectified images were generated with respect to the epoch of the ALS and GNSS surveys, irrespective of the overall planimetric displacements associated with the massive earthquakes, as well as the potential mountain uplift. This approach benefits from eliminating the possible systemic error (bias) due to tectonic displacements since 1956 by assuming a homogenous plate motion in the study area (i.e. roughly 13 km 2 , see Figure 1). Hence, a more consistent set of kinematic values for the studied landforms can be informed from a tectonically active region.

Kinematics Evolution and Recent Climatic Trends
Traditionally, rock glaciers have been considered less sensitive to climate change than glaciers (Barsch, 1996). This notion has been mainly held for the Andes, where long term subsurface temperature monitoring on mountain permafrost is lacking.
Yet, a few studies have observed permafrost degradation on Andean rock glaciers by geophysical investigations (Francou et al., 1999), short-term borehole measurements (Monnier and Kinnard, 2013) or by future climate warming scenarios (Drewes et al., 2018). Elsewhere in Europe, a plethora of studies have highlighted some symptoms of permafrost degradation, such as rock glacier acceleration Eriksen et al., 2018;Kellerer-Pirklbauer et al., 2018), destabilisation (Delaloye et al., 2013;Vivero and Lambiel, 2019;Marcer et al., 2021) or even collapse Marcer et al., 2020), mainly due to rising air temperatures and feedback mechanisms such as the availability of meltwater (Buchli et al., 2013;Cicoira et al., 2019). In our study area, despite their relatively small size, the kinematic results on the three most active landforms show velocity rates that are comparable to the rock glaciers in the northern Tien Shan region  and the Cordón del Plata range (Blöthe et al., 2020). The general acceleration trend since the 2000s that we observed with our results is in reasonable agreement with values reported elsewhere in the Northern Hemisphere (Hartl et al., 2016;Necsoiu et al., 2016;Eriksen et al., 2018;PERMOS et al., 2019;Kääb et al., 2021).
Regarding the atmospheric warming in this region (Carrasco et al., 2008;Falvey and Garreaud, 2009), our findings on the long-term kinematics evolution can be regarded as an indicator of environmental changes, in particular climatic and groundthermal conditions, but also the complex interaction between glaciers and mountain permafrost.
The short-term temperature data from the two shallow boreholes at Tapado 1 rock glacier (2 m depth, see Figure 1) suggests a slight cooling trend in the active layer during the 2010-2015 period (Figure 10), which is concomitant with the lower MAATs estimated at 4,500 m a.s.l during the same period and a short deceleration trend for the majority of the studied landforms ( Figure 7). This short cooling trend may be related to the prevalence between La Niña (dry and cold phase of the ENSO) and neutral ENSO conditions since 2009 (Garreaud et al., 2017). In addition, the active layer temperature variations revealed that the zero-curtain effect (Outcalt et al., 1990) was only present during the freezing May-June but not during the thawing conditions November-December. Commonly, the zero-curtain occurs during the onset of the freezing or thawing periods by delivering latent heat generated by water phase changes (French, 2017). Hence, the lack of a zero-curtain during thawing conditions can indicate low moisture content in the active layer (Zenklusen Mutter and Phillips, 2012). This could be explained by the premature depletion of the snow cover, occurring just before the temperature rises, and below the normal precipitation levels and high sublimation rates occurring during La Niña conditions (Ginot et al., 2006;Réveillet et al., 2020). Furthermore, the early removal of the thin snow cover facilitates the penetration of the winter cold in the active layer, potentially leading to a noticeable permafrost cooling. More research on the particular conditions of the Andean active layer and its relation to climate events is undoubtedly needed.

Evidence of Destabilisation Processes
Multitemporal analysis on the El Cachito rock glacier ( Figure 11) evidenced that the front toe of the rock glacier advanced about 20 m (∼1 m/yr) between 1978 and 2000, but only 6 m (∼0.3 m/yr) FIGURE 10 | Active layer temperature records (daily averages) at the two shallow boreholes (2 m depth) at Tapado 1 rock glacier (see Figure 1) Figure 11). During the first period, the lack of coherent surface topography precluded obtaining kinematic values just below the scarp, where the surface structure changed dramatically. Therefore, the surface velocities reported in Figure 7 could not capture the main destabilisation phase. The recent velocities (2014-2019, Figure 6) ranged from about 1.8 m/yr near the front (below the former scarp) to 0.35 m/yr near the upper central zone (above the former scarp). This form of reversed velocity profile has been observed on destabilised rock glaciers where extensive flow patterns are associated with convex slopes (Delaloye et al., 2013;Marcer et al., 2019) and might indicate that the destabilisation phase has not ended. Aside from the collapse of the Las Tórtolas rock glacier (20 km to the N, 4,500 m a.s.l, Bodin et al., 2012) or the sudden acceleration of the South Infiernillo Alto rock glacier after the artificial overload by waste deposits (330 km to the S, 3,900 m a.s.l, Valenzuela, 2004), there are no reports about destabilised or collapsed rock glaciers in the Andes. Nevertheless, a visual evaluation for the area surrounding the Tapado site (30°S) using Google Earth Pro revealed at least two other rock glaciers exhibiting destabilised morphologies (e.g. large transverse scarps and rapid terminus advance). The absence of detailed geomorphological maps and the lack of geophysical information hinders our ability to unravel the exact processes from these particular types of slope failures. Nevertheless, the destabilisation process in these rock glaciers stresses the need for more studies in the context of periglacial hazards in this part of the Andes (Iribarren Anacona et al., 2015;Vergara Dal Pont et al., 2020).

CONCLUSION
We have compiled one of the longest records of rock glacier kinematics in the Southern Hemisphere. The quantification of kinematics on different landforms was based on an assemblage of diverse remote sensing datasets covering the 1956-2019 period. Average rock glacier velocities ranged from 0.26 m/yr during the 1978-2000 up to 1.25 m/yr during the 2014-2019 periods. The recent rock glacier acceleration trend in this part of the Dry Andes mimics what has been reported in the Northern Hemisphere. UAV-derived datasets provided detailed short-term surface elevation changes and surface velocities at the northern frontal lobe of the Tapado complex, which were not observable with the other datasets.
Archival aerial photography, high-resolution satellite images and UAV surveys are suitable for determining the surface flow of different cryospheric components located in remote areas. The combination of UAV surveys and SfM photogrammetry can help document rapid cryospheric changes in a relatively low-cost approach. Furthermore, improved landform documentation can be achieved with the assessment of past and recent kinematics information. Also, our contribution can further support the work of the IPA Action Group rock glacier kinematics and inventories, which want to establish rock glacier kinematic as an Essential Climate Variable (ECV) in the framework of the Global Climate Observing System (GCOS) of the World Meteorological Organization (WMO). In this context, the appropriate analysis of archival aerial photography can help to reduce the paucity of long-term systematic measurements on rock glacier kinematics and fill the gap between Northern and Southern hemisphere permafrost observations.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
SV designed the study and wrote the manuscript under the supervision of CL, and SV performed the two UAV campaigns, processed the GNSS data, historical photographs and satellite imagery. DF-B, and XB helped with GNSS and archival aerial photography processing. SV, SM, NS, BR, and XB performed data acquisition during several field trips. All authors contributed to editing the manuscript.

FUNDING
SV recognizes the financial support of the Leading House for the Latin American Region (University of St. Gallen), grant number MOB1829.

ACKNOWLEDGMENTS
We greatly acknowledge the Geographical Military Institute (IGM) for the aerial photographs of Tapado and the Water Directorate of Chile (DGA) for the 2015 ALS data. We wish to thank the European Space Agency for the free provision of the Pléiades imagery through the Restrained Dataset project 41743. The original design of the surface rock glacier kinematics and the shallow borehole instruments, whose results are cited in the text, was initiated during the DGA-2010 funded programme 'Dynamics of rock glaciers in semiarid Chile' by Alexander Brenning, XB, and Guillermo Azócar (Pontificia Universidad Católica de Chile, Instituto de Geografía). DF acknowledges the Vicerrectoría de Investigación y Desarrollo, Universidad de Concepción, Postdoctorado VRID. We thank the editor, Matt Westoby, and two reviewers for their very useful comments on our manuscript.