We use the 2018 glacier outlines from Falaschi et al. (2019), drawn manually on ArcGIS on the basis of the 22 February 2018 Pléiades scene. Because the geodetic mass balance of rapidly receding glaciers is particularly sensitive to unaccounted area changes (Florentine et al., 2023), we revised the 2018 glacier outlines and adapted them by visual inspection of the 20 March 2023 scene.

For the purpose of feature tracking in optical imagery, we rely on rock glacier polygons available from Nansen and Belgrano catchment reports published in the National Glacier Inventory of Argentina (ING, Zalazar et al., 2020). We chose to use the ING outlines over the rock glacier inventory by Falaschi et al. (2015) because the ING follows a restricted mapping approach (Rock glacier inventories and kinematics, 2023). This allows us to emphasize the analysis on the rock glacier surfaces, dismissing rock glacier lateral margins and front talus, where image matching is less successful. To simplify the rock glacier name coding used in the ING in this study, we renamed the six investigated rock glaciers to Penitentes I-VI (Penitentes VI was not mapped in the ING), as per the Penitentes peak where the landforms are located.

We assessed surface velocities of 6 rock glaciers (accounting for a total area of 2.9 km²) for two epochs, i.e., between 29 March 2008 and 22 February 2018 (EP1), and between the latter and 20 March 2023 (EP2). For the initiation time of EP1 we used the 2.5 m resolution ALOS PRISM scene originally utilized by Falaschi et al. (2015) to produce a local rock glacier inventory. This ortho-image was generated using the DSM and Ortho-image Generation Software of the Japan Aerospace Exploration Agency with no terrain-surveyed ground control points (GCPs). The reader is referred to Falaschi et al. (2015) for full technical details and validation.

In addition to the ALOS PRISM dataset, we used two sub-meter Pléiades orthomosaics and Digital elevation Models (DEMs) at 0.5 m spatial resolution and 2 m ground-sampling distance, respectively, acquired on 22 February 2018 and 20 March 2023, and available from the Pléiades Glacier Observatory program (Berthier et al., 2024; LEGOS and OMP, 2024). Tested globally over mountainous terrain, the low uncertainty level (˂1 m, Berthier et al., 2014) of the very-high-resolution (VHR) Pléiades stereo imagery has allowed for robust assessments of glacier elevation changes at short (even annual to seasonal) time steps in various mountain ranges (e.g., Belart et al., 2017; Beraud et al., 2023; Berthier et al., 2024; Falaschi et al., 2023a).

The PGO standard procedure allowed us to assess glacier geodetic mass balance in Monte San Lorenzo between 22 February 2018 and 20 March 2023. The PGO workflow includes DEM coregistration, differencing and correction of systematic biases (see Berthier et al., 2024 and references therein). During further processing of the PGO data, we defined as outliers those cells in the elevation change grid maps exceeding ± 100 m and removed them accordingly. In DEM differencing, any data voids in the source DEMs are inevitably inherited by the derived elevation change grids. For our 2018–2023 Pléiades-derived elevation change grids in Monte San Lorenzo, data voids accounted for 1.5% of the glacier area. Within glacier ice area only (i.e., not over rock glaciers), we followed a global hypsometric approach (McNabb et al., 2019) to fill in these data voids, applying a fifth-degree polynomial fit to the elevation changes on 50 m elevation bins. Once all the above corrections were implemented, we derived time-averaged glacier- and region-wide elevation change rate values (dh/dt) at pixel scale.

The volumetric change Δv (m³) over the investigated time interval dt is then calculated as the product of the individual pixels dh _i (summed over the total pixel number covering the glacier area) and the cell size r of the dh/dt grid: ∆ v = r 2 × ∑ d h i

We then calculated the annual glacier-wide mass balance Δm (in m w.e. a⁻¹) using a density conversion factor ρ of 850 ± 60 kg m^-3. This value is technically sound in geodetic assessments considering periods longer than 5 yrs and significant volume changes (Huss, 2013): ∆ m = ∆ v × ρ Ᾱ × d t − 1 where Ᾱ is the mean glacier area calculated over the intervening time dt.

To calculate the volumetric mass balance uncertainty and overall uncertainty in mass budget, we follow the method described in Falaschi et al. (2023b), which uses an error propagation approach, considering the uncertainties in glacier area, elevation change and the volume to mass conversion factor. A full description of the calculation of the involved uncertainties and other consideration of biases can be found in Supplementary Section S1.

We derived surface velocities of 5 rock glaciers from 30 March 2008 to 22 February 2018 (EP1) and 6 rock glaciers from 22 February 2018 to 20 March 2023 (EP2) using the Environmental Motion Tracking software (EMT Version 0.93 - Schwalbe and Maas, 2017) and a feature tracking approach in optical imagery as described in Blöthe et al. (2025). While using the original resolution of the Pléiades imagery for tracking surface velocities for EP2, we resampled the ALOS PRISM and Pléiades imagery for EP1 to a common resolution of 2.5 m. For computational efficiency, we subdivided the study area into a northern and southern section, containing the rock glaciers Penitentes I to IV, and Penitentes V and VI, respectively (Supplementary Figure S1). To enhance rock glacier surface structures such as ridges and furrows for feature tracking in optical data recorded by different sensors, we applied a near anisotropic orientation filter during pre-processing (see Van Wyk de Vries and Wickert, 2021; Blöthe et al., 2025 for details). The workflow initially matches image patches from a reference to a secondary image utilizing cross-correlation, followed by a least squares matching (LSM) algorithm to achieve sub-pixel accuracies. For feature tracking, we defined object points as spaced 2.5 (EP2) and 5.0 m (EP1) in x and y directions within grids based on rock glacier outlines. We discarded object point displacements where correlation coefficients lay below 0.75. To account for positional mismatches between image pairs, EMT calculates an affine transformation model based on a set of 1,000 object points that were randomly distributed outside rock glacier polygons. To quantify the residual positional errors, we tracked a second set of 1,000 randomly distributed object points after applying the affine transformation. The Limit of Detection (LoD; Supplementary Table S2, Supplementary Figure S2) was then defined as the median of residual mismatches between image pairs, and is used to discern significant from non-significant displacements.

Recently, the Rock Glacier Inventory and Kinematics working group of the International Permafrost Association developed the Moving Area (MA) concept to identify areas of spatially homogeneous velocities within rock glaciers (Rock glacier inventories and kinematics, 2023). The MAs are used to assign a consistent kinematic attribute to a given rock glacier (active, transitional, relict), based on a predefined classification scheme. On the basis of the derived displacements, we transformed them to surface velocity rasters and subsequently classified them into the velocity classes as defined by the RGIK initiative (Rock glacier inventories and kinematics, 2023). Because the velocity field distribution of a single rock glacier can be highly variable, we set 0.01 km² as the MA minimum area, and considered rock glaciers as actively moving when at least one MA showed surface velocities above the corresponding Limit of Detection.

The region-wide geodetic mass balance (71 glaciers, 132.34 km²) around Monte San Lorenzo (Figures 1A, B) over the 2018–2023 period was −1.49 m ± 0.16 w.e. a⁻¹, and the overall volume loss amounted to 225.9 ± 17.8 × 10⁶ m³ a⁻¹ (Supplementary Table S2). To provide a sense of homogeneity in glacier mass balance, the non-area weighted average was −1.41 ± 0.16 w.e. a⁻¹, and the standard deviation (spread in mass balance among individual glaciers) was 0.73 w.e. a⁻¹. Among the larger sampled glaciers (i.e., ˃1 km²) Río Lácteo showed the most negative mass budget (−3.02 ± 0.22 w.e. a⁻¹). Significantly, the mean elevation change over 50 m elevation intervals (Figure 1D) shows thinning through the full altitudinal range of glacier ice in Monte San Lorenzo (450–3,700 m). Contemporarily, the total glacier area shrunk from 132.3 km² to 125.9 km² (∼1% a⁻¹).

Regarding rock glacier elevation changes, we found that on average, three sampled rock glaciers (Penitentes III-V) showed slightly positive elevation changes (up to +0.27 ± 0.88 m at Penitentes IV) between 2018 and 2023, whereas the remaining three (Penitentes I, II and VI) underwent surface lowering overall (Table 1), with a negative peak of −0.46 m at Penitentes VI.

Our feature tracking results indicate active movement exceeding the Limit of Detection for all investigated rock glaciers (Figure 2). For EP1, the Limit of Detection falls between 0.11 and 0.14 m a⁻¹, while being lower for EP2 with Limit of Detection values falling between 0.09 and 0.10 m a⁻¹. We identified 12 and 18 Moving Areas with velocities higher than the Limit of Detection for EP1 and EP2, respectively (Table 1), which account for 22.9% and 39.4% of the total rock glacier area. While no Moving Areas with velocities ˃1 m a⁻¹ were detected in EP1, Penitentes V attained the maximum velocity (1.1 m a⁻¹) during EP2, and included the (single) Moving Area with the highest velocity class (1–3 m a⁻¹), covering 14.1% of its total area. Velocities in the range of 10–30 cm a⁻¹ were the most abundant in absolute number (n = 6 and 11) and area (13.3% and 31% of the total rock glacier area) in both EP1 and EP2. Mean velocities of individual rock glaciers ranged from 0.20 m a⁻¹ to 0.47 m a⁻¹ in EP1 and from 0.14 m a⁻¹ to 0.52 m a⁻¹ in EP2, while median velocities varied from 0.16 m a⁻¹ to 0.42 m a⁻¹ in EP1 and from 0.14 m a⁻¹ to 0.43 m a⁻¹ in EP2.

From a sample of the 18 largest glaciers (covering 95% of the total glacier area), all but one show more negative mass budgets in 2018–2023 than during 2000–2018. Dedo Soto Glacier (0.9 km²) is the only exception, with a much more negative mass balance in 2012–2018 (−1.90 ± 0.16 m w.e. a⁻¹; Falaschi et al., 2019) than in 2018–2023 (−0.87 ± 0.13 m w.e. a⁻¹). The region-wide 2018–2023 mass budget (−1.49 ± 0.16 m w.e. a⁻¹) was twice as negative as in 2000–2012 (−0.69 m w.e. a⁻¹) and nearly three times more negative than in 2012–2018 (−0.58 ± 0.19 m w.e. a⁻¹). This reflects the close relation between summer air temperatures and glacier mass balance (Falaschi et al., 2019; Supplementary Figure S1 in this study). Furthermore, the 2018–2023 mass balance of the three largest valley glaciers was almost twice as negative as the 1958–2018 average (Falaschi et al., 2019; Supplementary Table S2). Although some elevation gains are observable in the upper parts of Río Oro, Río Lácteo and San Lorenzo Sur glaciers (Figure 1C), we put forward that these are most probably owed to DEM artefact in almost vertical terrain that are later transferred to the resulting elevation change grids.

Our 2018–2023 geodetic results provide the most up-to-date glacier mass change assessment at the periphery of the Patagonian Icefields, but confirming recent trends elsewhere in the Patagonian Andes remains challenging. The closest ongoing glaciological mass balance program, conducted by the ING since 2013, monitors Glaciar de los Tres in the Fitz Roy massif (∼200 km south of Monte San Lorenzo). There, unpublished data indicate a similar doubling of mass loss rates between 2018 and 2023 (−1.19 m w.e. a⁻¹) compared to −0.64 m w.e. a⁻¹ during 2013–2018 (WGMS, 2023; Pierre Pitte, personal comm.).

The analysis of temperature and precipitation records from the Balmaceda weather station (180 km South of Monte San Lorenzo, but representative of the overall climate conditions at Monte San Lorenzo, see Falaschi et al., 2019) reveals an annual temperature increase of 0.12 °C per decade and a precipitation decrease of 28 mm/d between 1963 and 2024 (Supplementary Figure S1). The mass balance of the San Lorenzo glaciers with the longest existing record (San Lorenzo Sur and Río Lácteo) appear nevertheless to be driven mainly by summer (Oct-Mar) temperatures (0.17°C/d), as shown by the Spearman Rank Order correlations in Supplementary Figure S1. While both winter and summer precipitation have decreased, the correlation between annual precipitation and glacier mass balance is not significant (r_s = 0.7 at p < 0.05). Because further increases in air temperatures are projected in the area around the Patagonian Icefields (Bravo et al., 2021), it is critical to actively monitor ongoing changes in the cryosphere in the region.

As outlined in Section 4.1, Penitentes III–V exhibited slight elevation gains, whereas, Penitentes I, II, and VI experienced slight elevation losses (Table 1). We stress here, however, that none of these changes are statistically significant (i.e., greater than the uncertainty). Nonetheless, it is telling to us that rock glaciers with an average positive elevation change are found in narrow gullies stripped of glacier ice, where rocky material can be avalanched from the surrounding slopes. Rock debris can detach from the headwall and linger on the surface of rock glaciers, later to be carried down the valley (Lehmann et al., 2022). This rock glacier configuration and associated nourishment mechanism are consistent with observations from other sites in the Andes, where positive elevation changes have been attributed to debris accumulation (Robson et al., 2022). On the contrary, rock glaciers with average surface lowering are comparatively larger and located in wider valleys with an upslope glacier connection. We put forward that during their interaction with glaciers, likely during the Little Ice Age, these rock glaciers might have accumulated significant amounts of ice in their upper sections. As ice was later covered by debris and is presently melting, this configuration led to negative elevation changes, similarly to what has been observed in other alpine settings (e.g., Vivero and Lambiel, 2024; Wee et al., 2024).

A closer look into the ‘raw’ (i.e., uncorrected for data voids) rock glacier elevation change maps, however, reveals that individual rock glaciers do show substantial changes over specific areas (Figures 3A–I). Unlike glaciers, which consistently thinned across all elevations (Figure 1C), some rock glaciers (Penitentes I, II, VI) experienced surface lowering (up to 10 m) only in their uppermost parts, where debris-free ice is present (Figure 3-a-). Below this, elevation changes remain within ±0.5 m, except in thermokarst ponds and ice-cliff depressions typical of debris-covered ice (Figure 3-b-). Further downslope, some elevation change maps reveal the downward movement of longitudinal and transverse furrows and ridges (Figure 3-e-). At the terminus, the advancing front of Penitentes I shows elevation gains exceeding 2 m (Figure 3-c-), and resembles the conveyor caterpillar-like advance mechanism of rock glaciers (Kääb and Reichmuth, 2005) probably unrelated to glaciers. Occasionally, rock glacier fronts display localized surface lowering (which we interpret as partially collapsed or eroded fronts), alongside adjacent areas of elevation gains, corresponding to the material mobilized and deposited immediately underneath (Figure 3-f-).

Ice-debris complexes are ice-debris landforms in close proximity to glaciers and may include several elements such as rock glaciers in various stages of activity, ice-cored moraines and also debris-covered stagnant ice bodies. Elevation changes in ice-debris complexes in regions exhibiting negative glacier mass balance often show significant surface lowering due to sub-debris ice melt in their upper parts and a greater stability towards their fronts (Bolch et al., 2019). The observed elevation change patterns in Monte San Lorenzo resemble those of ice-debris complexes with actively moving rock glacier-type fronts, as documented in the Andes and other mountain regions in the World (Monnier and Kinnard, 2017; Meng et al., 2023). ‘Typical’ rock glaciers, lacking clear connections to debris-covered or debris-free ice, rarely exhibit surface lowering of magnitudes of up to 20 m per decade (Groh and Blöthe, 2019; Cusicanqui et al., 2021; Fleischer et al., 2021, 2023). While such ice-debris complexes are well-documented in the drier central Andes (Blöthe et al., 2021; Villarroel et al., 2021; Robson et al., 2022), their presence in the more humid Central Patagonia had only been suggested based on geomorphological interpretation of high-resolution imagery (Falaschi et al., 2015; Zalazar et al., 2020). When observing landforms of analogous characteristics in the Jeinemeni district (∼80 km north of Monte San Lorenzo), Selley et al. (2019) attributed their potential origin to former glaciers that experienced highly negative mass balances since the Little Ice Age and are presently transitioning into rock glaciers (Monnier and Kinnard, 2015; 2017). While observed velocities and elevation changes support this, detailed geophysical analyses of internal structure and ice content are needed (Bolch et al., 2019; Petersen et al., 2020; Villarroel et al., 2021) to better understand the formation and evolution of these landforms in the Patagonian Andes. This could, in turn, provide insight into the geomorphological evolution of deglacierising mountain regions under future climate change.

From a methodological perspective, the high vertical precision of Pléiades-derived DEMs allows for detecting elevation and geomorphological changes in rock glaciers and ice-debris complexes, even when these changes are smaller than glacier thinning processes. While higher-resolution techniques like UAV surveys or Terrestrial Laser Scanning (Abermann et al., 2010; Zahs et al., 2019; Vivero et al., 2022) can capture finer rock glacier dynamics at seasonal or annual scales, very-high-resolution satellite photogrammetry remains effective for pluri-annual to sub-decadal investigations. Thus, Pléiades data are increasingly being used to monitor rock glaciers and ice-debris landforms in addition to glaciers, particularly in semi-arid regions where these landforms coexist (e.g., Robson et al., 2022; Cusicanqui et al., 2024; Wood et al., 2025).

By comparing the rock velocity maps from EP1 and EP2 (Figures 2C,F,I), we found that only a small proportion of the rock glacier area exhibited velocity changes that exceeded the EP1-EP2 difference Limit of Detection (Supplementary Table S2). Merely 4% of the total rock glacier area (and 10% of the area successfully tracked in both epochs) shows velocity increases ˃LoD. Similarly, velocity decreases ˂LoD were observed in less than 4% of the total rock glacier area and 9% of the common tracked areas. An exception is Penitentes V, where 21% of the total area and 29% of the successfully tracked areas showed velocity increases ˃LoD.

Despite significant warming documented in the Monte San Lorenzo area (Falaschi et al., 2019; Supplementary Figure S2), and a significant acceleration in glacier loss revealed in our data, rock glacier velocities remain stable since 2008 (within our detection limit), though our sample size remains somewhat small. While long-term data on rock glacier kinematics in the Andes is scarce and sample sizes vary across previous studies, findings remain mixed—some suggest widespread acceleration, while others indicate largely unchanged velocities (Vivero et al., 2021; Blöthe et al., 2025; Cusicanqui et al., 2024). The above findings are not surprising, since permafrost reacts in a much delayed and attenuated way and constitutes a longer-term indicator to climatic changes compared to glaciers (Kääb et al., 2007). Specifically, increasing air temperature, reduced ground-insulating snow cover, and concomitant increase in liquid water availability are regarded as the primary factors affecting permafrost temperatures and the deformation rate of frozen debris (Cicoira et al., 2021; Kääb and Røste, 2024). In addition, and unlike glaciers, which primarily respond to climate forcing, rock glacier dynamics is influenced by complex energy transfer processes and thermo-hydro-mechanical properties that affect heat conduction and storage in the ground (Wirz et al., 2016; Pruessner et al., 2018; Cicoira et al., 2019a). More so, thermal inertia and insulation effects of surface layers (e.g., snow) contribute to the delayed response of permafrost to climatic changes (Murton, 2021).