Investigating seasonal and multi-decadal water/ice storage changes in the Murtèl rock glacier using time-lapse gravimetry

Rock glaciers are important features of many alpine hydrological systems. Although their seasonal release of water enhances the resilience of alpine headwater catchments to climate change, measurement of their internal water and ice storage changes remains a challenge. Recent technological and methodological advances have enabled novel applications of time-lapse gravimetry (TLG) to estimate subsurface storage changes. Here, we present the first application of TLG on a rock glacier. We measure seasonal (July–September) changes in gravity at the Murtèl rock glacier (Upper Engadine, Switzerland). We employ drone-based photogrammetry to correct for surface mass changes in the form of snow. We also compare the Bouguer anomaly of our 2024 surveys with those from a pioneering 1991 gravimetry study. The seasonal results reveal spatial differences in active layer thaw, with estimated ice storage loss ranging between 11–64 cm water equivalent, while the multi-decadal results suggest zonal decreases in permafrost ice storage. Our study provides new insights into rock glacier–groundwater processes and illustrates how TLG can be employed to measure cryospheric and hydrogeological processes in permafrost and periglacial landforms.

1 Introduction

Our changing climate is imparting profound shifts to the hydrological regime of high-mountain regions (Hock et al., 2022). Temperature, evaporative losses, and precipitation variability are increasing, which could lead to more frequent droughts and reduced streamflow in dry years and thus threaten downstream water security (Schaffer et al., 2019; Hoelzle et al., 2020; Barandun et al., 2020; Arenson et al., 2022). As glaciers recede (Hugonnet et al., 2021) and precipitation in the form of snow decreases (Gottlieb and Mankin, 2024), the dynamics of many high-mountain regions are shifting from glacial to periglacial regimes (Haeberli et al., 2017). In this context of reduced surface snow and ice volumes and increasing risk of precipitation deficits, the hydrological buffering capacity of belowground components, namely ice-rich mountain permafrost and groundwater (Woo, 2011; Hayashi, 2020; Somers and McKenzie, 2020), is becoming increasingly important for water security and ecology (Millar et al., 2015; Reato et al., 2021). This is especially the case in water-stressed regions like semi-arid Central Asia (Karthe et al., 2015), the Himalayas (Jeelani et al., 2024), and the Dry Andes (Schaffer and MacDonell, 2022), but it is also becoming increasing relevant for other rapidly changing mountain regions such as the European Alps.

Intact rock glaciers (RGs), common in many high-mountain ranges, are ice-rich permafrost landforms that store and release water over various timescales ranging from seasonal to millennial (Jones et al., 2019). Their thick, seasonally thawed active layer (AL) insulates the ground and slows the melting of ground ice (Wakonigg, 1996; Humlum, 1997; Arenson et al., 2022; Amschwand et al., 2024). This renders permafrost ice more robust against climate change than glaciers over a long, multi-decadal timescale (Arenson et al., 2022; Haeberli et al., 2006). Three types of processes summarise the hydrological significance of RGs. First, seasonal storage and freeze–thaw of water/ice in the AL act as a seasonal buffer by building ground ice in autumn, winter, and spring (i.e., by trapping in-blown snow or refreezing snowmelt) and melting it during the thaw season. Second, permafrost ice in the permafrost body (beneath the AL) is slowly released due to permafrost degradation as intact (ice-bearing) RGs slowly transition towards a relict (ice-free) state. Since intact RGs are a common periglacial landform, the amount of belowground permafrost ice, as estimated from RG inventories and empirical area–volume scaling relations, is substantial (Rangecroft et al., 2015; Azócar and Brenning, 2010; Brenning and Azócar, 2010; Brenning, 2005; Jones et al., 2018b, a). Additionally, frozen talus slopes also contain ground ice (Mathys et al., 2024). In semi-arid, weakly glacierised catchments, water equivalent (w.e.) volumes stored in RGs can exceed glacier ice volumes or may do so in the future (Jones et al., 2018a; Bodin et al., 2010; Janke et al., 2017; Azócar and Brenning, 2010). Third, water is stored in the unfrozen fine-grained sediments of intact and relict RGs and interacts while flowing through or beneath RGs (storage–release, routing, and chemical alteration/mineralisation of water) (Azócar and Brenning, 2010; Corte, 1976; Burger et al., 1999; Jones et al., 2019). The storage capacity of liquid water increases with progressing RG degradation, where the available pore space increases at the expense of the lost ground ice (Wagner et al., 2016).

Ground-ice volumes in the AL are ostensibly much smaller compared to those of permafrost ice; however, they are the sites of rapid exchange between the atmosphere and hydrosphere. Marchenko et al. (2012) reported that observed seasonal accumulation and melt of 40–60 cm of ground ice stored a substantial amount (<30 %) of the snowpack in the Northern Tien Shan. Halla et al. (2021) found interannual ice storage changes on the Dos Lenguas rock glacier (Dry Andes of Argentina) of $-\mathrm{36}\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$ (25 %–80 % of the annual precipitation) and $+\mathrm{27}\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$ (17 %–55 %) for 2016–2017 and 2017–2018, respectively. On the Murtèl RG, where our present study is performed, seasonal ice storage changes in the AL were found to be 15–30 cm w.e. (10 %–32 % of the annual precipitation) based on point-scale belowground stake measurements and AL energy budgets (Amschwand et al., 2025a).

Distributed, non-invasive, and portable near-surface geophysical methods are well suited for mountain permafrost monitoring. Recent decades have seen the emergence of the subdomain of hydrogeophysics (Binley et al., 2015), wherein geophysical techniques are developed and applied to address hydrological and hydrogeological questions, and cryogeophysics (Hauck and Kneisel, 2008; Godio, A, 2019), to address cryosphere-related questions. As there is strong coupling of many hydrogeological and cryosphere processes (van Tiel et al., 2024), hydrogeophysical and cryogeophysical methods and applications have considerable overlap. Methods, such as electrical resistivity tomography (ERT), ground-penetrating radar (GPR), and active seismics, have found useful applications in alpine contexts. These include rockslide characterisation (e.g., Godio et al., 2006; Scapozza et al., 2011; Cody et al., 2020), aquifer delimitation (e.g., McClymont et al., 2010; Christensen et al., 2020; Halloran et al., 2023), glacier tomography (e.g., Egli et al., 2021; Dow et al., 2020; Church et al., 2019), and permafrost measurement (e.g., Maurer and Hauck, 2007; Draebing, 2016; Mewes et al., 2017). Geophysical methods can be considered either stationary or transient in their applications. Stationary measurements are suited to subsurface “mapping” and the differentiation of units or layers, while transient methods enable the measurement of spatio-temporal variations in water or ice content. The application of transient methods can reveal new, often quantitative, information on subsurface water and ice dynamics.

Gravimetry is a geophysical and geodetic (Jaramillo et al., 2024) set of techniques involving the measurement of g (acceleration due to gravity) across the Earth's surface. The GRACE and GRACE-FO missions (Rodell and Reager, 2023) measure transient g on a global scale; however, the >200 km spatial “footprint” of space-borne gravimetry is of little use for field-scale investigations. Stationary (steady-state) gravimetry at the land surface has a long history of providing information on geological interfaces. In permafrost studies, however, its use has been quite limited compared to other geophysical methods. Stationary gravimetry has been used to locate massive subsurface ice and ice-rich soils (Loranger et al., 2015) and to determine RG internal structure (Hausmann et al., 2007). At the Murtèl RG, a 1991 gravimetric survey was used to map the bedrock interface and spatial extent of layers observed in core samples (Vonder Mühll and Klingelé, 1994). Time-lapse gravimetry (TLG), sometimes referred to as “differential gravimetry” or “transient (micro)gravimetry”, involves measuring temporal changes in g at one or more locations. In alpine environments, use of the method is growing. McClymont et al. (2012) observed spatially variable Δg across a talus–moraine complex in the Lake O'Hara (Canada) headwater catchment during the snow-free period. Also over the snow-free period, Arnoux et al. (2020) recorded large decreases in g in the Vallon de Réchy (Switzerland) relative to a point located outside the catchment. Both of these studies used the Scintrex CG-5, which has since been superseded by the Scintrex CG-6, a significantly more stable and accurate relative gravimeter (Francis, 2021). At the Zugspitze (Germany), Voigt et al. (2021) used a continuous time series from a superconducting absolute gravimeter to investigate snowpack dynamics, while Gerlach et al. (2017) performed a similar investigation on the Vernagtferner glacier (Austria). TLG has not yet been explicitly used to measure permafrost or periglacial processes.

Here, we present the first TLG study of an active RG. At the well-documented Murtèl RG (e.g., Haeberli et al., 1998; Vonder Mühll et al., 2000; Hoelzle et al., 2002), we combine repeat measurements of gravity with repeat aerial imagery in order to determine the change in gravity attributable to ice melt in the AL (AL thaw) over the thaw season. Additionally, we calculate and compare the 2024 Bouguer anomaly (BA) data with those obtained in 1991 by Vonder Mühll and Klingelé (1994). We interpret these results in combination with point measurements (Amschwand et al., 2025b, a) to reveal spatial patterns in AL dynamics and multi-decadal permafrost loss. Finally, we discuss the potential of TLG for monitoring ice/water storage changes in the mountain cryosphere and give research perspectives on exploring permafrost–groundwater connectivity.

2 The Murtèl rock glacier

Murtèl RG (WGS 84: 46°25^′47^′′ N, 9°49^′15^′′ E; CH1903+/LV95: 2'783'080, 1'144'820; 2620–2700 m a.s.l. (metre above sea level); Fig. 1a) is located in a north-facing cirque in the Upper Engadine, a moderately continental, rain-shadowed high valley in the southeastern Swiss Alps (Fig. 1a). Mean annual air temperature (MAAT) is −1.7 °C, and mean annual precipitation (MAP) is ∼900 mm (Scherler et al., 2014). The RG covers an elevation range from 2620 (base of front) to 2720 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (transition to talus) (Fig. 1b–d). The talus slopes (2720–2800 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) connect the active RGs to the headwalls that rise up to 3165 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ The small (30 ha) non-glacierised Murtèl catchment consists of an unvegetated rock glacier, talus slope debris, and bedrock. Murtèl is located at the permafrost margin, and only its forefield (2600–2620 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) is permafrost free (Schneider et al., 2012, 2013). The rock glacier sits in a bowl-shaped, glacially overdeepened bedrock depression that has been mapped with a multi-geophysical approach combined with borehole logs (Vonder Mühll and Klingelé, 1994).

Figure 1 Locations of the gravity survey points and webcam photos from the survey dates. (a) Locations of the Maloja and Silvaplana points. Inset: location of the study site within Switzerland. (b) Locations of the local points on and in the vicinity of Murtèl RG. (c) Webcam image (courtesy of the Swiss Permafrost Monitoring Network PERMOS, University of Fribourg) of the RG during the first survey. (d) Webcam image of the RG during the second survey. Both webcam images have approximate locations of points MURTEL02 to -07 overlaid (MURTEL01 and -08 not visible). Background maps from Swisstopo.

The lobate Murtèl rock glacier is ca. 300 m long, 180 m wide, and covered by a coarse, blocky AL (debris mantle). Geophysical investigations revealed that the AL thickness varies according to the surface micro-topography, from ∼2 m in the furrows to ∼5 m beneath the ridges (Vonder Mühll and Klingelé, 1994; Vonder Mühll et al., 2000). Thus, the permafrost table shares the surface furrow-and-ridge micro-topography, albeit with attenuated relief. Fine material increases towards the AL base but is overall sparse. The permeable coarse, blocky AL appears to have a high permeability and a very low water retention capacity (Springman et al., 2012). The Murtèl permafrost body between the seasonally thawed coarse, blocky AL and bedrock (at ∼50 m) consists of three distinct layers (Vonder Mühll and Haeberli, 1990; Haeberli, 1990; Arenson et al., 2002): (1) massive ice, sparsely sand- and silt-bearing (3–28 m, supersaturated with over 90 % ice content); (2) a layer of ice-saturated frozen sand (28–32 m), accommodating ca. 60 % of the total/surface displacement (shear horizon); and (3) ice-saturated debris (32–50 m; 40 % ice). Three boreholes, all located within 30 m distance, intercepted this three-part stratigraphy but also revealed small-scale material differences laterally (e.g., variable ice/sand content, lenses) and thermal anomalies suggesting non-diffusive heat transfer and intra-permafrost water flow (Vonder Mühll, 1992; Arenson et al., 2010). The extent and ice content of the ice-rich permafrost body is well known (Vonder Mühll and Haeberli, 1990; Vonder Mühll, 1992; Vonder Mühll and Holub, 1992; Vonder Mühll and Klingelé, 1994; Haeberli and Vonder Mühll, 1996; Vonder Mühll et al., 2000). Electrical resistivity tomography (ERT) and seismic refraction tomography data (Hauck, 2013; Arenson et al., 2010; Mollaret et al., 2019) indicate that the ice-rich permafrost core has an extent of approximately 150×300 m², amounting to a water volume equivalent (WVEQ) of $\sim \mathrm{1.5}\times {\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{\mathrm{3}}$.

3 Methods: gravimetry

3.1 Basic principles of time-lapse gravimetry

Time-lapse gravimetry (TLG) is the measurement of acceleration due to gravity, g, at two or more points in time, usually at fixed locations, in order to determine Δg(t). More precisely, TLG measures changes in the vertical component of gravity, which is affected by changes in the distribution of mass:

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }g\left(t\right)=G\iiint \mathit{x}\cdot \widehat{z}{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }(\mathit{x},t)}{|\mathit{x}{|}^{\mathrm{3}}}}{\mathrm{d}}^{\mathrm{3}}x,\end{array}$$

where G is the universal gravitational constant, $\mathrm{6.674}\times {\mathrm{10}}^{-\mathrm{11}}\mathrm{N}{\mathrm{m}}^{\mathrm{2}}{\mathrm{kg}}^{-\mathrm{2}}$; $\widehat{z}$ is the vertical unit vector; and Δρ is the change in density at location x relative to the measurement point. Δg measurements are generally reported in units of mGal (${\mathrm{10}}^{-\mathrm{5}}\mathrm{m}{\mathrm{s}}^{-\mathrm{2}}$) or µGal (${\mathrm{10}}^{-\mathrm{8}}\mathrm{m}{\mathrm{s}}^{-\mathrm{2}}\approx {\mathrm{10}}^{-\mathrm{9}}\mathit{\text{g}}$). Individual surveys using a portable relative gravimeter generally include multiple points and, through repeat measurements at one or more reference stations, form one or more closed loops. After corrections to the raw data (Sect. 3.3), the residual effect on g due to mass distribution changes is revealed. Local mass loss due to ice melt or decreased groundwater storage in the vicinity below a measurement location will result in a decrease in gravity between two repeat surveys. As different mass change processes can result in the same Δg values at a given location, interpretation must be informed by other information (Sect. 6.3).

3.2 Time-lapse gravimetry surveys

We carried out two microgravimetry surveys on 8 July and 11 September 2024 (Δt=65 d) using a Scintrex CG-6 Autograv (Scintrex, 2018). A 60 s integration time was used, and, to ensure stability, the device was set to measure until three consecutive measurements within a tolerance of ∼2 µGal were recorded. The average of this stable 180 s duration measurement period was used in the below-discussed analysis. To enable time-lapse measurements and to correct for instrument drift, both surveys started and terminated at absolute gravimetric reference point no. 1018 (Swisstopo, 2024), a first-order point ($g=\mathrm{980}\mathrm{224}\mathrm{831}\pm \mathrm{3}\mathrm{\mu }\mathrm{Gal}$) located on a bedrock outcrop at the Maloja Pass. In addition to the MALOJA reference point, nine points were surveyed during both surveys: Silvaplana church, two locations (MURTEL01 and -08) adjacent to the RG, and six locations on the RG (MURTEL02 to -07) (Fig. 1). The number of points and their locations were chosen so as to maximise coverage across the RG and to limit the survey duration, ensuring a better drift correction. The total duration of each survey was under 5 h 20 min. Future surveys at Murtèl – or, indeed, at other RGs – could include more measurement locations, requiring multiple survey loops and a longer survey duration. We marked survey points with paint and used a Leica GNSS DGPS¹ to survey all points to < 2 cm absolute accuracy. All survey points were located directly on bedrock or boulders. One challenge with gravimetric surveys on RGs is finding surfaces that are sufficiently horizontal so that the gravimeter can be levelled. Additionally, the height of the gravimeter above the ground surface was measured with a measuring tape to millimetre precision at each location.

3.3 Gravity corrections

Several standard corrections were applied to the gravimetric data: tilt, internal temperature, height difference, Earth tides, and drift. The first two, tilt and temperature, are automatically applied by the gravimeter to each individual measurement. Tilt was maintained in the calibrated <10 arcsecond range during all measurements; thus, the internal corrections were ≲1 µGal. The gravimeter was in a stable internal temperature state during the surveys, having being supplied with continuous power for several days prior to the surveys, thus ensuring internal temperature corrections were within the calibrated range. Corrections for height, Earth tides, and drift were all applied manually. During repeat measurements at a given location, the exact height of the device can vary slightly due, in part, to offsets caused by levelling. Gravity decreases with elevation (free-air gravity gradient):

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial g}{\partial r}}=-{\displaystyle \frac{\mathrm{2}G{M}_{\mathrm{E}}}{r^{\mathrm{3}}}}\simeq -\mathrm{3.1}\mathrm{\mu }\mathrm{Gal}{\mathrm{cm}}^{-\mathrm{1}},\end{array}$$

where M_E is the mass of the Earth and r is the distance to its centre of mass. While the true value of the local vertical gravity gradient (VGG) varies over Earth's surface, especially in complex terrain, 3.1 µGal cm⁻¹ is adequate for this correction of a few mm.

Earth tides, whose magnitude can reach ∼300 µGal, are periodic variations in g primarily due to the relative positions of the Sun and Moon with respect to the Earth. Through tidal harmonic databases, Earth tides can be calculated to a level of precision and accuracy that far exceeds gravimetric measurement accuracy. Modern gravimeters implement these corrections internally; however, to ensure the highest possible accuracy, we chose to manually apply these corrections using the current state-of-the art Kudryavtsev (2004) database, which contains over 28 000 wave groups. Through a Python script, we used the pygtide package (Rau et al., 2022) to calculate the Earth tides for each individual measurement to sub-nGal precision. These corrections, which ranged from −62.5 to +81.7 µGal for the July survey and from −45.9 to +41.1 µGal for the September survey, were applied to the measured data. The resulting datasets were used for the calculation of, and correction for, instrumental drift, which occurs in all relative gravimeters. The Scintrex CG-6 has a low instrumental drift of $<\mathrm{200}\mathrm{\mu }\mathrm{Gal}{\mathrm{d}}^{-\mathrm{1}}$ and $<\mathrm{20}\mathrm{\mu }\mathrm{Gal}{\mathrm{d}}^{-\mathrm{1}}$ residual drift (Scintrex, 2018), with independent observations of $<\mathrm{1.5}\mathrm{\mu }\mathrm{Gal}{\mathrm{h}}^{-\mathrm{1}}$ (Francis, 2021). By comparing the repeat measurements at MALOJA (Fig. 1a), we calculated and corrected for total drift. Any residual effects of transient atmospheric pressure (Klügel and Wziontek, 2009) are also indirectly corrected for during the drift calculation. Finally, the corrected data were referenced to the MALOJA reference point to enable calculation of the temporal difference, Δg, in the two surveys.

3.4 Bouguer anomaly (BA) calculations for comparison with 1991 gravimetric data

The gravimetric investigations of Vonder Mühll and Klingelé (1994), carried out in Summer 1991, provided information on the internal structure of the Murtèl RG. Raw data, such as the gravity station coordinates and non-terrain corrected gravity values, from the 1991 survey are unfortunately no longer available. Thus, while terrain corrections are not necessary for time-lapse gravimetry studies, since the data in the Vonder Mühll and Klingelé (1994) study are Bouguer anomalies (i.e., terrain corrected), we also perform these corrections to our gravimetric data to enable comparison. Furthermore, while our gravimetric stations MURTEL05, -06, -03, -07, and -08 correspond to “Profile 3” (Vonder Mühll and Klingelé, 1994; Vonder Mühll, 1993), we did not return to the exact locations of the 1991 stations.

Calculation of Bouguer anomalies (g_BA) involves four corrections:

$$\begin{array}{}\text{(3)}& g_{\text{BA}}=g_{\text{measured}}-[{\mathit{\gamma }}_{\mathrm{0}}+c_{\text{FA}}+c_{\text{BP}}+c_{\text{terrain}}],\end{array}$$

where γ₀ is “normal” gravity as calculated from the ellipsoid, while c_FA, c_BP, and c_terrain are the free air, Bouguer plate, or terrain corrections, respectively. While methods for all of these corrections (also sometimes referred to as “reductions”) are extensively documented (e.g., Nowell, 1999; LaFehr, 1991), there are several variations in the possible formulas and methods for each of these corrections. The methods reported by Vonder Mühll and Klingelé (1994), and also detailed in Vonder Mühll (1993), included an erroneous form of the 1967 international gravity formula in a second-order expansion (Caputo and Pieri, 1968), a linear $(-\mathrm{0.3086}\mathrm{mGal}{\mathrm{m}}^{-\mathrm{1}})z$ free-air correction, and a linear Bouguer plate correction of $(\mathrm{4.19}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{mGal}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{kg}}^{-\mathrm{1}})\mathit{\rho }z$ (where $\mathit{\rho }=\mathrm{2670}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$). Their reported terrain corrections included two inner Nagy (1966) prism zones and two outer line-of-mass zones to 167 km. Our application of these methods did not result in BAs in their reported range from −157.0 to −155.4 mGal. Such issues related to reproducing historic BAs published in the Swiss Gravimetric Atlas (Swisstopo, 1994) have been encountered and documented by others (Bandou et al., 2022). Swisstopo, the Swiss Federal Office of Topography, has also confirmed these issues related to incomplete documentation of the processing of historic datasets. We thus apply modern standard corrections in order to enable comparisons of BAs relative to others in the same survey. Specifically, we use the full GRS80 Somigliana formula (Moritz, 1980) for normal gravity, the Wenzel free-air correction (Wenzel, 1996) with ETRS89 latitudes, and the Talwani Bouguer plate formula (Talwani, 1973) with a density of 2670 kg m⁻³. Terrain corrections were performed using the Swisstopo GRAVNIV software with the publicly available 0.5 m Swiss digital terrain model swissALTI3D, as discussed in Bandou et al. (2022) and further detailed in Appendix 3A of Bandou (2023).

3.5 Digitisation of the 1991 survey data

As discussed in Sect. 3.4, the only sources of the 1991 survey data currently available are the graphical data in Vonder Mühll and Klingelé (1994) and Vonder Mühll (1993) (both of which contain the same information). We used a digitisation tool to extract, directly from the figures, the positions of the 1991 gravimetry survey points and the values of the calculated BAs and relative position along “Profile 3” data. Accuracy in the BA data digitisation procedure can be considered $\pm \sim \mathrm{10}\mathrm{\mu }\mathrm{Gal}$ and ∼2 m along the profile. Planar positions can be considered accurate to $\pm \sim \mathrm{3}\mathrm{m}$, although we note some discrepancy between planar positions on the original map (Vonder Mühll, 1993) and the distances along the original “Profile 3”. In the absence of other information, we present the 1991 survey data as originally published.

4 Methods: snow spatial distribution and corrections using aerial imagery

In order to distinguish between surface and subsurface cryosphere processes in our gravimetric analyses, we need to account for the effect of snow present during the first survey. To accomplish this, we used photogrammetry with UAV imagery.

4.1 Aerial imagery data acquisition

On both survey dates, we obtained aerial images using a compact DJI Mini Pro 4 with a 48 megapixel camera. The drone was piloted over the lower part of the RG to ensure coverage across the target zone. Due to the presence of overhead aerial cables and an operating gondola lift, extended zones were not surveyed. In total, 375 images were acquired during the July 2024 survey and 268 images during the September 2024 survey. As non-RTK² UAVs have poor positional accuracy, we also established and measured ground control points (GCPs) in order to improve outputs from the photogrammetry processes. We used a Leica GNSS RTK surveying system to measure the 3-D positions of eight GCPs, colocated at the gravimetric survey points MURTEL01 to -08 (Fig. 1), to 0.1 cm precision and <2 cm accuracy.

4.2 Photogrammetry

Digital surface models (DSMs) and orthomosaics were created for both survey dates. We used the photogrammetry software Pix4Dmapper (v4.9.0, 2023, Pix4D S.A., Switzerland). The average ground-sampling distance was 2.04 cm for the July survey and 3.45 cm for the September survey. A high number of median key points per image (>72 000 for both surveys) were determined, each with a minimum of four matches across multiple images. We manually georeferenced six GCPs (MURTEL02 to -07) that could be located in multiple images with high precision and accuracy. The georeferenced point clouds were used to produce DSMs and orthomosaics with 5 cm resolution. We used noise filtering and inverse distance weighting, but no surface smoothing, for the DSM processing. Comparison of the final results revealed only very small offsets at the sub-pixel scale.

4.3 Snow density measurements

In order to calculate the effect of snow on measured gravity, snow density needed to be measured. We measured the density using a portable digital balance (precision 0.1 g) and a coring tube (inner diameter 5 cm). The volume and mass of 18 samples from different locations at various depths up to ∼1 m were measured. The measured sample sizes ranged from 137–736 cm³ and from 76–410 g. Measured densities ranged from 504–587 kg m⁻³, a typical range for névé (Fierz et al., 2009). The average snow density was determined to be $\mathrm{550}\pm \mathrm{27}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$.

4.4 Calculation of the snow gravity effect

The gravitational effect of the snow, present in the July survey and absent in the September one, at each of the survey points was calculated through a multistep process: (1) creation of a snow raster mask based on the July orthomosaic, (2) calculation of the height difference between the DSMs, (3) fusion of these data to create a 3-D differential snow volume model, and (4) forward calculation of the vertical gravity effect at the gravity survey points. Python, implementing standard packages (numpy, matplotlib, etc.) and rasterio, was used for these processes. Throughout, manually clipped DSMs and orthomosaics were used, limiting the analysis to snow-covered regions with high-quality DSM data for both surveys.

We first created a spatial mask for the snow-covered pixels using the July 2024 orthomosaic (Sect. 4.2). This georeferenced image was converted to a luminance (L) greyscale format as per the Rec. 709 formula:

$$\begin{array}{}\text{(4)}& L=\mathrm{0.2125}R+\mathrm{0.7154}G+\mathrm{0.0721}B,\end{array}$$

where R, G, and B are the red, green, and blue pixel intensities, respectively. The luminance image was blurred using a uniform, centred filter of dimensions 1.2 m×1.2 m (25 pixels×25 pixels). This was necessary in order to avoid isolated high-luminosity pixels being identified as snow. A simple 55 % threshold was applied to the processed image, with the exception of the low-albedo snow in the southeast corner of the surveyed zone, where a 38 % threshold was used. This process resulted in satisfactory identification of snow-covered regions, while avoiding misidentification of non-snow-covered regions such as high-luminosity boulders. The height difference between the July and September DSMs was evaluated by a simple difference between the two. The above-discussed snow mask, as well as a 5 cm thickness lower bound, was applied to this raster data. The final result was a 2-D map of snow thickness (Fig. 2), which, when combined with either the July or September DSM, provides 3-D information on the change in snow volume distribution.

Figure 2 Snow depth and the digital surface models (DSM) calculated via photogrammetry and manual fixed point referencing. (a) The depth of the snow present during the July 2024 survey and absent in the September 2024 survey. Inset (b) is a histogram of the area covered by the snow with depths in 10 cm intervals. (c) The DSM for the 8 July 2024 survey. (d) The DSM for the 11 September 2024 survey. Only areas with sufficient coverage by both surveys are shown. The black dots indicate the locations of our gravimetric survey points (Fig. 1).

To calculate the decrease in gravity due to the melting snow, we performed forward gravity modelling using a modified version of Gravi4GW (Halloran, 2022) that implements the Nagy prism (Nagy, 1966) rather than equivalent point masses. A correction for the height of the proof mass (6.58 cm above the bottom of the gravimeter, Scintrex, 2018) relative to our surveyed coordinates (Fig. 1b) was also applied, although the effect was minor for all points (0.0–0.3 µGal).

Additionally, when comparing the digital surface models from the two surveys, we discovered a new large boulder, ostensibly originating from the unstable zone located adjacent and above (southwest) the RG. We estimated the boulder's mass to be $<\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{kg}$. The nearest gravimetry point, MURTEL04, was located ∼25 m away and ∼1 m lower than the boulder's centre of mass. The straight-line gravitational pull of this boulder on MURTEL04 was thus ≲1 µGal, and the vertical component (i.e., that which we measure in time-lapse gravimetry) was <0.1 µGal. Thus, we conclude that this new rock mass can be ignored in our analyses.

5 Results

5.1 Seasonal gravity variations

In order to correct for the gravity effect of the snow present in the July 2024 survey, the net gravitational effect on the survey points MURTEL01 to -08 was calculated (Sect. 4.4). This vertical Δg effect of the snow ranged from −0.58 (MURTEL04) to −11.39 µGal (MURTEL06) (Table 1). Only the vertical component of snow-induced gravity change has an effect on measured gravity (Eq. 1 and Sect. 6). While a significant amount of snow was present (Fig. 2), the majority of its gravitational effect is non-vertical. For reference, 1 m³ of snow (with a density of 500 kg m⁻³) at a horizontal distance of 10 m, centred 1 m below a given point, would have a vertical gravity effect of $\sim \mathrm{3.6}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{\mu }\mathrm{Gal}$. We also note that while snow-covered regions extending hundreds of metres from the RG could not be surveyed, these zones are generally at higher elevations than the survey points. The gravitational effect of the snow in those lateral zones would thus be small and positive (McClymont et al., 2012). The limited spatial coverage of the DSMs would thus have the effect of underestimating net mass loss in the AL and subsurface water. Finally, the statistical uncertainty in the measured snow density was 4.9 %, which translates to an uncertainty ranging from ±0.03 (MURTEL04) to ±0.56 µGal (MURTEL06).

Table 1 Δg results after processing (Sect. 3.3) and Δg−Δg_snow after accounting for the effect of surface snow mass Δg_snow.

* Incomplete DSM coverage of surrounding terrain.

Accounting for the snow, net gravity changes were negative at all eight survey points (Table 1), indicating a local decrease in subsurface mass. With the exception of MURTEL04, located near the western flank of the RG, all net gravity changes were in the approximate range of −5 to −15 µGal. Along the central profile, Δg decreases in magnitude from −15.2 µGal near the lower rooting zone (MURTEL05) to −5.2 µGal near the front. MURTEL08, located on bedrock ∼20 m from the RG front, experienced a decrease in gravity of −10.6 µGal. Results from the two lateral survey locations, MURTEL02 near the east flank (−10.3 µGal) and MURTEL04 to the west (−29.2 µGal), differ significantly, although limited DSM coverage may have a minor influence on the snow-corrected value at MURTEL04. Overall, while there is a spatial trend along the central line, the results point to significant heterogeneity in the subsurface processes. For first-order interpretation, the Bouguer plate approximation (BPA), an infinite flat plane where a decrease of 2.38 cm in water level equivalent equates to a −1 µGal change, provides a rough starting point, implying decreases in subsurface water storage from 12–70 cm of water equivalent (see Sect. 6.1.1 for an in-depth discussion and quantitative interpretation).

For uncertainty estimation, we follow the same logic as previous TLG investigations (McClymont et al., 2012; Arnoux et al., 2020), noting that these two studies used Scintrex CG-5 devices, which have significantly worse temperature and tilt stability and greater drift rates compared to the CG-6 device used in our study (Francis, 2021). In terms of uncertainty, we assume three significant and uncorrelated sources of uncertainty, in addition to the before-discussed uncertainty in the snow and rockfall correction (Table 2): precision, drift, and elevation. Precision groups together uncertainty related to vibrations, tilt corrections, and internal temperature corrections (Sect. 3.3). As determined from standard deviation of the three repeats whose average was used to determine the final values, we observe values ranging from ±0.2 to ±2.4 µGal with an average of ±0.9 µGal for the July survey and ±0.05 to ±1.1 µGal with an average of ±0.6 µGal for the September survey. Adding these uncertainties together in quadrature for each survey point, we obtain precision uncertainties in Δg of ±0.7 to ±2.5 µGal. We assume an uncertainty associated with the impact of non-linear drift of ±4.4 µGal based on survey duration and device specifications. Finally, uncertainty in Δg associated with uncertainty in the height of the device is assumed to be ±1 µGal as we manually measured the height of the device to ≲0.3 cm precision (Eq. 2) in addition to surveying. Calculations of these uncertainties, combined with the snow and rockfall correction uncertainty for individual points, range between ±4.6 to ±5.2 µGal and thus we arrive at a total estimated uncertainty in the snow cover-corrected Δg of ±5 µGal.

Table 2 Δg uncertainty budget.

* The total uncertainty is calculated from the square root of the sum of the squares of the individual uncertainty components.

5.2 Comparison of the 1991 and 2024 gravity surveys

As discussed (Sects. 3.4 and 3.5), even though most TLG applications do not require normal gravity, free air, Bouguer plate, and terrain corrections, we needed to calculate and apply these to calculate the Bouguer anomaly for comparison with the digitised 1991 survey data. Furthermore, due to incomplete documentation of the processing steps in Vonder Mühll and Klingelé (1994) and lack of the original DSM files, absolute differences between the published 1991 values of our BA values differ by an offset which cannot be known exactly (Sect. 3.4). While uncertainty of the 1991 survey data is not extensively discussed in Vonder Mühll and Klingelé (1994) or Vonder Mühll (1993), the techniques and technology of the time clearly limited the achievable accuracy of the data. The LaCoste and Romberg model G gravimeter used in 1991 had a repeatability of 100 µGal and accuracy of <40 µGal. The higher density of points in that study (1 every ∼25 m in the central line profile) allowed for the calculation of a smoothed trend line which we use for comparisons to the 2024 data on a relative basis (Fig. 3).

Figure 3 (a) Bouguer anomalies (BAs) and (b) approximate relative positions. The smoothed trend line (Sect. 3.5) from the 1991 survey (Vonder Mühll and Klingelé, 1994; Vonder Mühll, 1993) and those calculated from our 2024 surveys (Sect. 3.4) are shown. Distance refers to that of “Profile 3” in the 1991 survey, starting from its uppermost point. Approximate distances of the 2024 survey points in proximity to the 1991 profile are shown. An additional point near “Profile 1” in the 1991 survey is not shown: −152.447 and −152.476 mGal in July/September 2024 vs. $\sim -\mathrm{156.15}\mathrm{mGal}$ in 1991 (Vonder Mühll and Klingelé, 1994). Note that, due to lack of documentation and differences in methodologies (Sect. 3.4), the values cannot be compared directly, although spatial trends and the relative values within each survey can.

By displaying the calculated BAs with a fixed offset of 3.9 mGal (Fig. 3), we can approximately match our off-RG forefront survey point, MURTEL08 (Fig. 1), to the 1991 smoothed line in order to facilitate visual comparisons of the relative values. Our Bouguer anomalies display a larger difference between the forefront region (MURTEL08) and front and central RG zones (MURTEL07, -03, and -06) than in the 1991 survey. In contrast, the difference between the forefront and lower rooting zone at the transition to the talus slope (MURTEL05) is similar to that observed in the 1991 survey. Comparing with the smoothed historic data, the relative difference between the forefront and the front/central zones appears to be approximately 110–190 µGal greater in 2024 than in 1991. Similarly, comparing MURTEL08 to MURTEL04 (Fig. 3), the difference between the BAs at these points increases by ∼200 µGal between the 1991 and 2024 surveys. Our other survey points, MURTEL01 and -02, do not overlap with any of the 1991 survey points. Their BA values are −152.241 µGal (July) and −152.253 mGal (September) and −152.550 (July) and −152.564 mGal (September), respectively. Comparisons of data from our surveys and the Vonder Mühll and Klingelé (1994) survey must be interpreted with caution due to the unquantified uncertainties in the 1991 survey data and the differences in the locations of our survey points. As such, we stress that the 1991–2024 comparison provides qualitative, relative insights rather than precise, quantitative measurement of Δg.

6 Discussion

6.1 Seasonal processes and gravity

6.1.1 Seasonal storage changes in the active layer

The corrected July–September Δg values (Table 1 and Fig. 4) result from the integrated effect of mass distribution changes surrounding the survey point, weighted by distance $r^{-\mathrm{2}}\cdot \widehat{z}$ (Eq. 1). Corrections for change of mass on the surface, in the form of snow (Sect. 4), allow for local subsurface processes to be isolated. The net negative gravity changes Δg indicate a decrease in total water storage (ground ice and liquid water) – that is, net export of ice/water out of the rock glacier and underlying aquifers. While uncertainties exist (see discussion below), we attribute the dominant Δg signal over the thaw season to ice storage changes in the AL. However, we caution that, on their own, the observed seasonal decreases in g can not be attributed exclusively to active layer thaw or to groundwater storage decreases. Our reasoning is based on previous observations, specifically belowground stake measurements of the moving ground-ice table and ground-ice melt inferred from calculated AL energy budgets (Amschwand et al., 2025a). In order to estimate the water equivalent (w.e.) of this change in ice mass, we employ the Python tool GRAVI4GW (Halloran, 2022), which takes into account the non-horizontal and non-planar topology of a changing unconfined groundwater table or, in our case, that of the ground-ice table (AL base). The sensitivity of gravity changes, Δg, to total water storage changes,

$$\begin{array}{}\text{(5)}& {\mathit{\beta }}_z(\mathit{x},h_{\text{eff}})={\displaystyle \frac{\partial \mathit{g}}{\partial h_{\text{eff}}}}\cdot \widehat{z},\end{array}$$

varies with the effective depth of the groundwater/ice table h_eff and is influenced by the surrounding topography (Halloran, 2022). β_z varies with terrain slope and curvature just as the local vertical gravity gradient (VGG) does (Zahorec et al., 2024) and is highest on convex topographic features (mounds, ridges). A related concept is topographic admittance (e.g., Chaffaut et al., 2022; Voigt et al., 2021), which is most relevant for thin surface water accumulation (e.g., in the form of snow). The influence of small-scale topographic features decreases, while the effective radius of influence or “footprint” increases, with increasing h_eff (Halloran, 2022). Based on field observations (Amschwand et al., 2025a), PERMOS borehole temperature series (Noetzli and Pellet, 2023), and geophysical investigations (Vonder Mühll and Klingelé, 1994; Vonder Mühll et al., 2000; Mollaret et al., 2020), the AL depth can be assumed to be in the range of 2–5 m. Using this range of values for h_eff, the measured Δg values (with the ±5 µGal uncertainty) translate to an estimated water equivalent of AL ice melt Δh_we between surveys in the range of 15–40 $\mathrm{cm}\mathrm{w}.\mathrm{e}.$ for all on-RG locations except MURTEL04 (Figs. 4 and 5).

Figure 4 The 2024 seasonal gravity and total water storage changes at survey points MURTEL02 to -08. Total water storage (in $\mathrm{cm}\mathrm{w}.\mathrm{e}.$) calculated with h_eff=4 m (uncertainty ranges shown in Fig. 5).

Figure 5 Calculated seasonal (8 July–11 September 2024) belowground total water storage change Δh_we for a range of hypothetical h_eff values (see Eq. 5) and measured Δg values (thick lines) ±5 µGal (dotted lines) for the on-RG survey points MURTEL02 to -07 (see Fig. 4 for locations of survey points). The off-RG survey points, MURTEL01 and -08, are included for completeness.

TLG measurements are fundamentally ambiguous because the gravitational effect depends on the magnitude, distance, and relative position of the mass change relative to the survey point (Eq. 1). Two features of TLG are noteworthy (and discussed further in Sect. 6.3): gravity does not discriminate between water in its solid or liquid state, and it is in places weakly sensitive to the depth where the water storage changes occur (Chaffaut et al., 2022; Halloran, 2022). As a consequence, different groundwater and ice storage changes – whether localised and shallow or extended and deep (or any combinations thereof) – can produce the same Δg signal that reflects the total water storage change below ground (i.e., excluding snow). This is shown in Fig. 5, where total water storage change Δh_we is plotted as a function of depth h_eff to the (hypothetical) water/ice table. While other subsurface processes (Sect. 6.1.2) may play some role, the inferred Δh_we of 15–40 $\mathrm{cm}\mathrm{w}.\mathrm{e}.$ compares well with the results of Amschwand et al. (2025a) (15–30 $\mathrm{cm}\mathrm{w}.\mathrm{e}.$), suggesting that the snow-corrected on-RG Δg values are dominantly a signal of ice storage changes in the AL. This interpretation is most convincing at the mid-frontal survey points MURTEL06, MURTEL03, and MURTEL07 (Fig. 4), which are in the vicinity of direct AL thickness measurements (Amschwand et al., 2025a). Due to its sparse fine-material content, the Murtèl AL has a small retention capacity for liquid water, and supra-permafrost groundwater export is negligible on a seasonal scale. This interpretation is supported by the flashy outflow during the thaw season. After precipitation events in late summer, the supra-permafrost aquifer drains rapidly within ∼5 d, and the rock glacier springs dry out (Amschwand et al., 2025a). Furthermore, electrical resistivities of the Murtèl active layer are high, and injecting current has been notoriously difficult to measure, also suggesting a low water content (Noetzli and Pellet, 2023). Murtèl is a coarse, blocky rock glacier, and the effect of changing supra-permafrost groundwater levels is potentially important on rock glaciers, having more fine material in the AL. Under the assumption that the dominant mass redistribution effect stems from AL ice melt, the uncertainty in AL storage change water equivalent, Δh_we, is in the range of 10–14 $\mathrm{cm}\mathrm{w}.\mathrm{e}.$ for survey points on the RG. This uncertainty stems from both the uncertainty of the final snow-corrected Δg values and the uncertainty in the local AL depth below the surface. The sensitivity–depth relation β_z (Eq. 5) differs between survey points and varies with terrain slope and curvature. The Δh_we uncertainty stemming from h_eff (as opposed to from Δg measurements) is influenced by topography and can be calculated beforehand to optimise the survey point locations. Survey points on moderately convex features (ridges and plateaus) generally give higher sensitivity to storage changes, β_z, that varies little with effective depth h_eff. This is illustrated by the mid-frontal survey points MURTEL06, MURTEL03, and MURTEL07, where β_z (and consequently h_we) varies little with h_eff. There, the Δh_we uncertainty stems mainly from Δg uncertainty (Table 2) rather than that of h_eff. With a β_z (Eq. 5) weakly sensitive to depth h_eff, total storage change estimates h_we are only weakly impacted by the uncertainty in AL depth or the possible occurrence of deeper thaw processes.

6.1.2 Other seasonal-scale processes

In principle, Δg is sensitive to more than just ice loss by AL thaw. Other changes in total water storage, in both solid and liquid states, and mass movements can also affect gravity at the RG surface (Fig. 6). These include effects from snow, intra-permafrost melt and subsidence, groundwater export, and rock movements. Snowpack can significantly contribute to the Δg signal and must be accounted for by calculating and subtracting its differential gravitational effect (Sect. 4). Based on the drone-derived differential DSMs and the snow density measurements, the snow cover effect could be precisely calculated (within ±0.6 µGal, Table 2) and amounted to $\sim \mathrm{20}-\mathrm{50}\mathit{\%}$ of the measured Δg signal (Table 3), except at survey point MURTEL04. There, near the western rock glacier margin, the snow's vertical gravity effect was counterintuitively weak despite large snow depths locally exceeding 4 m in the nearby depressions (Fig. 2). The small gravitational effect of the snow occurred here, because mass losses above and below the survey point have opposing effects on the measured vertical Δg and have, effectively, partially cancelled each other out (melt at higher elevation increases g, while melt at lower elevation decreases g). However, limited extended DSM coverage towards the adjacent Marmugnun rock glacier precludes a conclusive uncertainty analysis for the correction at this point.

Figure 6 Conceptual longitudinal cross-section model (not to scale) displaying water/ice storage changes that occur in rock glaciers. While the structure and bowl-shaped base of Murtèl are integrated here, the displayed processes may occur in all rock glaciers. On the seasonal scale of our July–September 2024 TLG measurements, superficial melt of snow (which we account for in Sect. 4.4) and near-surface AL thaw (discussed in Sect. 6.1.1) are the dominant processes. On multi-year timescales, all of the processes (labelled in red), as well as rock mass movements (Sect. 6.2), may play a role. These processes may alter the value of g as measured at the surface.

Belowground water/ice storage changes (in excess of the pore space) can result in subsidence or heave. At Murtèl, however, subsidence (loss of excess ice) is $\le \mathrm{1}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$, an order of magnitude smaller than the seasonal ice storage change (Amschwand et al., 2025a). While the presence of liquid water in the permafrost core of Murtèl (and other rock glaciers) is well documented from borehole drilling (Vonder Mühll and Haeberli, 1990; Arenson et al., 2002; Bast et al., 2024), little is known about temporal changes in intra-permafrost water content. Based on the weak heat fluxes (Amschwand et al., 2025b) and the low hydraulic permeability of the ice-rich permafrost body (Springman et al., 2012), we consider effects of intra-permafrost thaw (that would decrease g) and changes in water content negligible on seasonal timescale compared to the substantial ice turnover in the AL.

Assuming low hydraulic conductivity and porosity in the crystalline bedrock (granodiorite), groundwater export from the base of the Murtèl RG, lying in an over-deepened bowl-shaped depression (Barsch and Hell, 1975; Vonder Mühll and Klingelé, 1994), could reasonably be expected to be rather limited on the seasonal timescale. Nonetheless, weathered bedrock at shallow depths may have significantly higher conductivities and porosities than at depth, allowing for seasonal storage changes that may affect g to a measurable level. The decreases in g at the off-RG points suggest that groundwater storage has decreased. While both MURTEL01 and MURTEL08 are likely to be influenced by RG active layer thaw, MURTEL08 in particular (Fig. 5h) suggests that groundwater storage decreases have also occurred. There is, however, no borehole intersecting the bedrock; thus, it is impossible to draw firm conclusions about groundwater storage changes beneath the Murtèl RG, where known seasonal AL ice thaw occurs and has been directly measured. Based on the observation that the RG springs dry out in hot and dry summer periods when melt rates are highest, Amschwand et al. (2025a) found that a significant proportion of Murtèl meltwater does not leave the catchment as surface streamflow and thus most likely recharges the sub-permafrost aquifer. Halla et al. (2021) reached the same conclusion for a rock glacier in the Dry Andes. The influence of groundwater recharge at the RG base is nonetheless impossible to verify without bedrock-intersecting piezometric levels, hydraulic characterisation of the fractured bedrock, or detailed water balance measurements and calculations. This cryosphere–groundwater connectivity, whose relevance may be significant on longer timescales, is further discussed in Sect. 6.4.

As for non-water/ice mass redistribution phenomena, accumulation of rock mass above and near the survey points could effectively reduce g. Due to the inverse-distance-squared weighting (Eq. 1), mass movement effects on Δg would be most pronounced at survey points closest to the active root zone (talus slope) and to the nearby Marmugnun RG (MURTEL04), although this depends on the individual rockfall trajectories. The gravitational pull of a large rockfall boulder deposited between the two time-lapse survey dates has been calculated from the two DSMs (Sect. 4.4), and these processes are deemed insignificant at the seasonal scale, although this depends on the exact locations of survey points and rockfall trajectories. In this study, we do not account for processes occurring in neighbouring regions, although we note that decreasing mass at elevations above our survey points would be positive and thus have the effect of “masking” the change in gravity imparted by AL ice melt and other subsurface ice and water mass losses in the RG.

6.2 Multi-decadal processes and gravity

Comparison of the 1991 survey data (Vonder Mühll and Klingelé, 1994) and our measurements must be nuanced by the technological and methodological factors discussed in Sect. 5.2. Nevertheless, some interpretations and hypotheses can be set forth. While over a single thaw season, or part thereof, melt of ground ice in the AL is the dominant mass redistribution process, over multi-year timescales – in this case 33 years – there are several processes that may significantly affect g.

The starkest change in the spatial trends of the BA (Fig. 3) is the relative decrease in gravity in the central and frontal regions (MURTEL06, -03, and -07) as compared to both the forefront measurement point (MURTEL08) and that in the lower rooting zone (MURTEL05). Although MURTEL08 is located ∼25 m from the 1991 survey profile, the others are all located <5 m from it. Thus, even if we ignore the forefront point, a clear difference in trend can be observed where the BA is now lower relative to that of the rooting zone than it was in 1991. Mountain permafrost regions are highly dynamic environments characterised by mass-moving processes (rock/debris) and total water (water, snow, ice) storage changes acting on different magnitudes, rates/timescales, and depths below ground, i.e., on the surface, in the AL (supra-permafrost), in the permafrost body (intra-permafrost), and beneath (sub-permafrost). The Murtèl periglacial cirque is typical in this respect (Müller et al., 2014). Here, we discuss processes related to both rock mass movements and water/ice storage changes that could be relevant for decadal changes in gravity in alpine periglacial settings (Table 3).

Landform and rock mass movements involve gradual or episodic movement of mass. Active rock glaciers undergo creep at rates of a few cm yr⁻¹ to m yr⁻¹ (Delaloye et al., 2010; Noetzli and Pellet, 2023; Manchado et al., 2024; Louis et al., 2025), and destabilised rock glaciers can move at rates in the $\sim \mathrm{cm}{\mathrm{d}}^{-\mathrm{1}}$ range (Marcer et al., 2021; Hartl et al., 2023). Even at slowly deforming RGs, such as Murtèl, the effects of creep on mass distribution, and thus potentially g, will be significant on multi-year timescales. Furthermore, RG creep will eventually render it impossible to return to the same survey points for repeat surveys as the point in space may end up being within the rock mass or high above the land surface. Although Murtèl exhibits comparatively slow creep rates of $\lesssim \mathrm{10}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$, amounting to $\sim \mathrm{2}-\mathrm{4}\mathrm{m}$ of horizontal displacement over 33 years, the displacement alone is very unlikely to account for the 1991–2024 differences. We note that shifting these 2024 BA points by $\sim -\mathrm{30}\mathrm{m}$ along the x axis in Fig. 3 would result in an overlap with the 1991 trend line. Rockfall and other landslide-type events also redistribute mass. The interplay of the effects of rock debris from headwalls above rooting zones and those of RG creep is highly complex and would require careful and extended topographic surveying (i.e., with lidar or photogrammetry) at the time of each gravimetric survey. When conceptualising the gravitational effects of these, and, indeed, all mass redistribution processes, it is important to remember the $r^{-\mathrm{2}}\cdot \widehat{z}$ weighting (Eq. 1) of mass changes. All things being unchanged, if mass is removed from an elevation above the point and added at an elevation below a point, then gravity will increase. However, the effect of falling rock mass is dependent on the relative distance and vertical orientation relative to a measurement location, both before and after the rockslide episode. Finally, tectonic uplift may be important to consider for multi-decadal TLG surveys, although its effect would be uniform across such a small study area at the Murtèl RG. In the Bernina Alps, tectonic uplift rates are about 1–2 mm yr⁻¹ (Sternai et al., 2019); thus, tectonic uplift contributes no more than $-\mathrm{0.6}\mathrm{\mu }\mathrm{Gal}{\mathrm{yr}}^{-\mathrm{1}}$ (Eq. 2) or −20 µGal to the BA between 1991 and 2024.

Water storage changes, in both liquid and solid form, impart both seasonal (periodic) and decadal changes to g. Surface processes, such as snow accumulation and melt, cause periodic changes and are much more straightforward to measure than subsurface ones, although gravimetry has been demonstrated to be useful in its monitoring on short and long timescales (Breili and Pettersen, 2009; Voigt et al., 2021). Annual freeze–thaw cycles in the AL are also primarily periodic in their effects (Sect. 6.1.1). Decadal progressive deepening (or thinning) of the AL could potentially be measurable with TLG. In periglacial environments in general, permafrost degradation will have measurable effects on g, as will groundwater storage changes (see Sect. 6.4). In RGs, specifically, intra-permafrost melt of excess ice results in a net loss of mass and may result in an increase in bulk density. At a given point in space above the RG, these processes would decrease g, although the effect on BA would be more complex to discern due to the potential for compaction and also due to dynamic effects from compressive thickening or extensive thinning. Overall, our BA results suggest that, compared to 1991, there is now a relatively lower mass density in the central and frontal regions compared to the rooting zone. It is highly likely that both ice/water storage (Fig. 6) and rock movement processes have influenced these changes in the BA. Additionally, the periodic drying of frontal springs during melt periods (Sect. 6.1.2) supports the hypothesis that deep groundwater export may have also occurred. Future high-precision gravimetric measurements that can be compared in absolute terms to our 2024 surveys will enable the evaluation of long-term water/ice storage change processes within the RG.

Table 3 Types of mass distribution changes in mountain permafrost environments that potentially affect the surface value of g and thus are relevant for time-lapse gravimetry. Periodic processes (e.g., AL freeze–thaw) will dominate on short (seasonal) timescales, while other processes may have significant effects over longer timescales.

6.3 Potential of TLG for monitoring of permafrost and periglacial landforms

TLG offers clear opportunities for monitoring the various seasonal and long-term water/ice storage processes that occur in permafrost and periglacial environments. While gravimetry has existed for decades, modern portable devices, such as the Scintrex CG-6 used in our study, offer levels of stability and accuracy that provide new opportunities for monitoring subsurface processes whose effects may only be on the few µGal level. Gravity measurements are non-invasive, integrative, and distributed, and they are possible essentially anywhere that is accessible by foot. Application of TLG to determine absolute gravity changes requires a known reference point, which may not always be feasible. Switzerland has a dense network of maintained absolute gravimetric reference points (Swisstopo, 2024), and the Murtèl RG has the particularity of being rapidly accessible via cable car for much of the year. This advantageous situation is not always the case for alpine sites. Nevertheless, there are also multi-step strategies that can be employed to establish temporary local references, although they may require large travel distances and long times, depending on the location. These constraints can effectively prevent the carrying out of absolute TLG surveys. In these cases, relative TLG, which involves using a reference point at which the value of g is not known absolutely and may vary with time, may be an option. Relative TLG, wherein a fixed Δg offset between surveys is not known, still provides data on changes in g at survey points relative to one another and has already provided useful knowledge in the growing field of alpine hydrogeology (Arnoux et al., 2020).

The non-discriminatory sensitivity of gravimetry to mass changes can be an advantage in permafrost or periglacial environments, where water co-exists in solid and liquid state, as it yields measurements of the total water storage change including ice and water. Without the need for petrophysical relations, forward modelling can be relatively straightforward. This makes gravimetry an ideal complement to other geophysical measurements that are selective to specific physical properties (e.g., electromagnetic or elastic properties). Conversely, attributing the total water storage change to individual supra-, intra-, and sub-permafrost ice/water storage changes (Table 3) is ambiguous and requires independent observations or a priori knowledge of the stratigraphy and hydraulic properties. Gravimetric surveys also require precise positional measurements of survey points to enable accurate temporal comparisons. For example, in the case of subsidence caused by excess ice melt, mass is lost beneath the survey point, which should decrease gravity. Simultaneously, the subsiding point moves closer to Earth's centre of mass, increasing gravity. Only over multi-year timescales can the two competing effects be accounted for, and the loss of the excess ice can only be quantified if elevation and VGG are known to a high precision. As many complex processes occur in periglacial landforms over seasonal, yearly, and longer timescales, TLG surveys at multiple timescales could provide valuable information to help constrain hydrological, structural, and rheological models of these geomorphological structures.

6.4 Permafrost–groundwater connectivity

As recently underlined in multiple recent reviews (Somers and McKenzie, 2020; Hayashi, 2020; van Tiel et al., 2024), alpine hydrogeological processes and, in particular, cryosphere–groundwater connectivity are insufficiently understood. In many alpine headwater catchments, RGs are a vital component of annual hydrological and hydrochemical dynamics (Brighenti et al., 2019). Permafrost ice, whether in RGs or other landforms such as frozen talus slopes (Mathys et al., 2024), is, in essence, immobile groundwater whose presence constitutes a water store and also modifies the subsurface hydraulic properties. Ice build-up and melt are, in turn, controlled as much by the ground thermal regime and the heat fluxes (Amschwand et al., 2024) as by water availability from precipitation or snowmelt, be it the seasonal turnover in the AL or permafrost thaw in the underlying rock glacier core. Semi-impervious permafrost bodies can act as aquitards (Pavoni et al., 2023), leading to rapid supra-permafrost runoff of intense precipitation and a “flashy” hydrograph typical of many permafrost-underlain catchments (Rogger et al., 2017). At Murtèl RG, however, the before-discussed hydrological measurements, stake measurements, and AL energy budgets suggest that up to 20–30 % of the snowpack is retained as AL ice and melts slowly enough to percolate through the semi-impermeable permafrost body. This “leaky” behaviour and the hypothesis that coarse, blocky permafrost landforms contribute to groundwater recharge by seasonally storing and routing winter precipitation to a sub-permafrost aquifer have been investigated at multiple sites (Amschwand et al., 2025a; Corte, 1987; Halla et al., 2021; Marchenko et al., 2024). TLG, as discussed in Sect. 6.3, is sensitive to all mass changes; therefore, the method would need to be combined with other techniques in order to resolve and differentiate storage changes in permafrost and groundwater.

Due to its extensive documentation and overdeepened bowl-shaped bedrock base, Murtèl would be an ideal site for further studies investigating permafrost–groundwater connectivity. The concave underlying bedrock profile may facilitate the infiltration of basal flow to the fractured bedrock (Fig. 6). Characterisation of the fracture network in the granodiorite bedrock would allow for a more comprehensive understanding of the cryosphere–groundwater processes occurring at Murtèl and their implications for the wider catchment. A multi-method approach, combining hydrological and geophysical measurements with numerical hydrological modelling would help us to understand the timings, volumes, and flow paths of water originating from AL ice melt and permafrost degradation.

7 Conclusion and outlook

Climate change is inducing profound changes in the mountain cryosphere and is rapidly modifying the hydrology of high-mountain watersheds. Water percolation across permafrost bodies, as well as intermittent storage in the form of ice, controls groundwater recharge in permafrost-underlain watersheds. These cryo-hydrogeological processes contribute to headwater baseflow, which is critical for aquatic habitats and water supply during droughts (Hayashi, 2020). Despite the importance of these phenomena, the complex links between the cryosphere and groundwater are poorly understood (van Tiel et al., 2024). There is thus a strong need for quantitative and spatially distributed methods for monitoring intra-annual, inter-annual, and long-term changes in subsurface ice content and water storage in these environments. Recent technical advances in hydrogeodesy (Jaramillo et al., 2024) have opened new avenues for investigating the connectivity of groundwater with the rapidly changing cryosphere. These advances complement established cryo-geophysical techniques (Hauck and Kneisel, 2008). Our application of time-lapse gravimetry (TLG) to estimate seasonal ice storage changes in the active layer of the Murtèl RG is the first use of this geophysical technique in a periglacial environment and has demonstrated the ability of the technique to spatially resolve AL thaw. Furthermore, despite the limitations of the historical data, our comparison with Vonder Mühll and Klingelé (1994) has shown the potential of multi-decade TLG measurements for permafrost degradation studies. With respect to our combined TLG–photogrammetry methodology and its application to an active RG, we observe the following:

1. TLG is an ideal method for spatially distributed quantification of mass movements and water/ice storage changes on the seasonal scale due to its non-invasiveness, portability, and spatially integrative sensitivity to mass changes.

2. The potential for long-term permafrost monitoring using TLG is promising, although comparison with historic gravimetry data is challenging unless the methodology is documented in detail and unless Δg changes exceed the potentially large uncertainties of historic data. The initiation of high-precision gravimetry–photogrammetry (or gravimetry–lidar) surveys at new sites would provide an invaluable baseline for the monitoring of permafrost degradation and its associated hydrological changes.

3. In an environment as dynamic as the mountain cryosphere, digital surface models (DSMs) from rapid photogrammetric surveys are valuable – in many cases necessary – for corrections of changing surface snow, névé, or ice coverage and rockfall deposits between surveys. They could also serve to correct for landform movements in long-term TLG studies or those rapidly deforming ice-rich permafrost landforms.

Through continued interdisciplinary collaboration across the geophysics, geodesy, geomorphology, cryosphere, and hydrogeology research communities, we anticipate that methodological innovations and diverse perspectives will result in a deeper, more quantitative understanding of total water storage change processes in alpine environments. In mountain permafrost and periglacial environments, TLG holds significant potential for impact on its own and would provide unique, complementary data in joint deployments alongside other transient geophysical, geochemical, and hydrological monitoring techniques.