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, orthomosaic and registration) and quality control steps were performed by Digimapas. After filtering, the final point density was 13.67 pts/m² with a vertical accuracy (RMSE) within 0.086 m (DIGIMAPAS, 2015).

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.

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 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, 2020). 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.

Two stereo panchromatic images from the GeoEye-1 and two tri-stereos Pléiades satellites were available for this study (Table 1). GeoEye-1 and Pléiades images were provided in level 1B and 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.

The UAV surveys were planned to acquire recent high-resolution images of the northern frontal lobe of the Tapado complex (Figure 1) and performed simultaneously with the GNSS surveys in January 2018 and March 2019. The 2018 and 2019 UAV surveys were conducted using Phantom 3 Pro and Phantom 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².

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).

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 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 the nearest meteorological station (La Laguna, 3,160 m a.s.l.) using the Chilean Climate Explorer (http://explorador.cr2.cl). In this study, we derived the MAAT from 1976 (beginning of the temperature measurements) to 2019.

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.

Figure 5 displays surface velocities of the Tapado site derived from a pair of Pléiades orthoimages between 2014 and 2019. This figure highlights the complex morphology and heterogeneous kinematic patterns of the Tapado complex. It also reveals the noticeably coherent kinematic behaviour of the Las Tolas and El Cachito landforms. Areas with steep slopes, snow cover and significant textural changes between 2014 and 2019 provided low-quality kinematic data and therefore were excluded from our analysis. However, surface kinematics on the debris-covered glacier revealed rather chaotic values due to drastic 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).

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 the most active rock glaciers is revealed during the same period.

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 high-resolution 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.

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.

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², see Figure 1). Hence, a more consistent set of kinematic values for the studied landforms can be informed from a tectonically active region.

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 (Kääb et al., 2007; 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 (Bodin et al., 2017; 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 (Kääb et al., 2021) 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 ground-thermal 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.

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) between 2000 and 2019. Note that during the same periods, the front toe of the Tapado 1 rock glacier advanced only 4 m (∼0.2 m/yr) and 5 m (∼0.25 m/yr), respectively. Also, a transverse scarp approximately 120 m long at an elevation of 4,390 m asl can be observed at the El Cachito rock glacier from the 1978 orthoimage (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).