Alpine rockwall erosion patterns follow elevation-dependent climate trajectories

Mountainous topography reflects an interplay between tectonic uplift, crustal strength, and climate-conditioned erosion cycles. During glaciations, glacial erosion increases bedrock relief, whereas during interglacials relief is lowered by rockwall erosion. Here, we show that paraglacial, frost cracking and permafrost processes jointly drive postglacial rockwall erosion in our research area. Field observations and modelling experiments demonstrate that all three processes are strongly conditioned by elevation. Our findings on catchment scale provide a potential multi-process explanation for the increase of rockwall erosion rates with elevation across the European Alps. As alpine basins warm during deglaciation, changing intensities and elevation-dependent interactions between periglacial and paraglacial processes result in elevational shifts in rockwall erosion patterns. Future climate warming will shift the intensity and elevation distribution of these processes, resulting in overall lower erosion rates across the Alps, but with more intensified erosion at the highest topography most sensitive to climate change. Rockwall erosion from glacier retreat, permafrost and frost cracking increases in intensity at higher elevations in a pattern that is sensitive to climate change, according to field observations, models and a rockwall erosion compilation in the European Alps.

T he topography of glaciated mountains reflects an interplay between tectonic uplift, crustal strength, and climate 1,2 . During glaciations, the glacial buzzsaw 3,4 can limit summit elevations, but glacial erosion mostly increases relief by carving deep valleys 5,6 . During and after deglaciation, non-glacial hillslope erosion processes, that tend to limit relief, takeover. Glacier retreat typically exposes steep, unsupported rockwalls that erode via paraglacial slope failure 4,7,8 . Paraglacial slope failures are directly conditioned by glacial activity [9][10][11] , and are prepared and triggered by glacial debuttressing [12][13][14] , internal stress redistribution 15,16 and seismic activity 17,18 following deglaciation. Periglacial processes operate on the same rockwalls, and include frost cracking and permafrost dynamics 3,4,7,8 . Frost cracking is a climate-driven 19 mechanical weathering process that breaks down rock 20 . Permafrost is a state of frozen ground, which can stabilise rock by increasing shear resistance 21 , but can also amplify frost cracking processes 22 and cryogenic stresses, which can, along with permafrost degradation, weaken rock masses 23 . These periglacial processes prime and trigger rockfalls 24 , contributing to rockwall erosion, and have been conceptualised as a relief-limiting frost buzzsaw 25 .
As paraglacial and periglacial processes cohabit similar topography, postglacial erosion of most alpine rockwalls is likely to be driven by a combination 4,7,8 of glacier retreat [26][27][28][29] , frost cracking [30][31][32] and permafrost processes [33][34][35] . Climate warming will modify the intensity and overlapping distribution of all three processes, and so anticipating trajectories of alpine topography evolution and rockfall hazard requires holistic quantification of these drivers. Investigation of these processes has tended to be in isolation and seldom at a landscape-scale, such that the relative contributions of these processes and their spatial distribution has hitherto never been quantified.
Here we provide a systematic analysis of postglacial rockwall erosion drivers by the integration of field and modelling experiments. We use empirically-validated models to evaluate the contributions of permafrost, frost cracking, and glacier retreat to rockwall erosion within an alpine basin. Our data indicate that paraglacial and periglacial processes combine to control climatedriven rockwall erosion in our research area and that both are elevation dependent, contributing to increased erosion rates at higher elevations. Our findings provide a potential explanation of rockwall erosion rate pattern that we compiled for the European Alps that is characterised by increasing erosion with increasing elevations. This relationship is amplified at catchments with mean rockwall elevations above 2500 m (all elevations refer to m above sea level). As these processes are linked to climate, climate warming will reduce periglacial and paraglacial erosion at lower elevations but intensify it at higher elevations. As a consequence, we propose that periglacial and paraglacial processes jointly exert a limit on postglacial summit heights, extending the existing frost buzzsaw hypothesis 25 .
At locations with mean rockwall elevations below 2500 m, the long-term erosion rates range from 0:252 þ0:148 À0:142 mm a −1 to 0:5 þ0:5 À0:15 mm a − 1 40,51,52,54,55 and are higher than the short-term erosion rates of 0:05 þ0:03 À0:02 to 0:225 þ0:075 À0:075 mm a −1 40,41 from equivalent elevations. An explanation for this may be that longterm erosion rates reflect a period of more intense periglacial processes and paraglacial adjustment 59 (consistent with conceptual paraglacial exhaustion curves) that would have existed immediately following cool Last Glacial Maximum (LGM) conditions 60 . Long-term erosion rates may also integrate a period of elevated post-glacial seismic activity, suggested to have increased mass movements in some locations 18 . In contrast, the short-term rates below 2500 m represent areas that are presently glacier-and permafrost-free and with relatively mild frost cracking 32 .
For rockwall elevations above 2500 m, which are potentially affected by permafrost (Fig. 1c) 43,45,46 or increased frost weathering within bergschrunds 47 , but such hypotheses have not been robustly tested and have tended to ignore paraglacial effects.
The time-scale dependant contrast in erosion rates above and below 2500 m suggests that elevation, though its influence on local climate, is a fundamental control on erosion processes. The patterns we observe are unlikely due to data preservation bias or the Sadler effect 61 (like has previously been observed in glaciated landscapes with short term rates exceeding long-term rates, interpreted as bias 62 ). The data we compiled reveal no consistent pattern of short-term erosion rates being greater than long-term rates (Fig. 1d), with instead long-term rates exceeding short-term rates at lower elevations. Sedimentary incompleteness can introduce bias 61 , but the long-term rockwall erosion rate data are largely from coarse-grained sedimentary archives such as talus slopes 40,[49][50][51][52][53][54]58 or rock glaciers 48,57 that are unlikely to have major sediment depletion gaps. There is limited opportunity for reworking and erosion of such coarse-grained deposits, which are located mostly within low-connectivity systems, such as disconnected glacial hanging valleys, and are net stores of postglacial sediment 63 . Furthermore, as our data compilation is from a single orogen we can exclude bias arising from comparing sites that experience markedly different tectonic regimes. In the following section, we evaluate the hypothesis that elevation is a first-order control on rockwall erosion rate patterns through its influence on periglacial and paraglacial processes.
Rockfall processes and potential influencing factors. To explore the relationship between elevation and rockwall erosion rate for the European Alps, we quantified and modelled periglacial and paraglacial processes and compared these to patterns of contemporary and post-LGM rockfall activity in an alpine valley. We chose the Hungerli Valley (Fig. 2) Supplementary Fig. 1). The rockwalls have been glacially eroded and experienced phases of glacial advance and retreat since the LGM. With periglacial conditions and a glacier terminus above 2850 m, the Hungerli Valley is representative of many alpine valleys in the European Alps (Fig. 1c). The rockwalls studied occur on the north and north-east facing slopes of Hungerlihorli peak and Rothorn cirque, respectively, with elevations from 2500 to 3277 m.
We detected 263 rockfall events in the valley between 2016 and 2019 using annually repeated terrestrial laser scanning. Rockfall volumes ranged from 0.1 to 159.4 m 3 (Fig. 3a), with a median of 0.26 m 3 , and 96% being smaller than 10 m 3 (Fig. 4e). The rockfall magnitude-frequency data exhibit a typical inverse relationship (N 2016-19 = 17.59 V −0.68 ; Fig. 3a), with a power law exponent, within the range (−0.48 to −0.77) of those previously published for alpine rockwalls 40,43,47,66 , that differs little between the three years (  Fig. 4a). The median rockfall volume scales positively with increasing elevation below the 3000-3100 m elevation band and negatively above this band (Fig. 4e), rather than following the distribution of available rockwall area proportionally (Fig. 4a). Most high-magnitude (> 10 m 3 ) Fig. 1 Rockwall erosion rates along glacier and permafrost distribution in the European Alps. a Compiled short-term, mid-term and long-term erosion rates in the European Alps along with the extent of LGM glaciation 144 , current glaciers 36 and permafrost distribution 38,114 (DEM source: European Digital Elevation Model). b Relationship between erosion rates and elevation. The linear relationship increases from r 2 =0.46 for all erosion rates to r 2 =0.63 (p < 0.001) when limited to catchment areas with mean elevation >2500 m. c Inset histogram: frequency of very small glaciers (<0.5 km 2 ; light blue), small glaciers (0.5 to 5 km 2 , blue), medium-sized glaciers (5 to 10 km 2 , turquoise) and large glaciers (>10 km 2 , dark blue) in the Alps. Main graph: Elevation of the terminus 36 of different sized glaciers (using same colours as for inset histogram) and permafrost distribution (purple line) 38,114 . d Comparison of short-, midand long-term erosion rates, to evaluate influence of integration period (i.e. time-scale) on calculated erosion rates. Error bars indicate the minimum and maximum range of rockwall elevation and erosion rates. For detailed erosion rates including used methods, see Supplementary Table 1. events occurred between 2900 and 3200 m, which encompasses most of the elevation range of the present-day glacier (Fig. 2).
As a preliminary exploration of rockfall patterns in the Hungerli Valley, we used a principal component analysis (PCA) to examine statistical relationships between modelled or inferred periglacial and paraglacial conditions, and topography. For each rockfall we related the deglaciation age (i.e. a paraglacial proxy; from historical aerial photos, moraine positions, and ice surface reconstructions), mean annual rock surface temperature (MARST; from a rock surface temperature driven permafrost model 67 ), elevation (i.e. a frost cracking proxy), slope angle (i.e. shear stress proxy), and the rockfall volume. The results show that the first two PCA dimensions explain 74.0% of the data variation (Fig. 3b), with a strong relationship found for deglaciation age (Eigenvalue 0.84) and MARST (0.72), and an inverse relationship with elevation (-0.85) on dimension 1 (PC1). The results suggest that all three elevation-dependent processes (deglaciation history, permafrost, and frost cracking) correlate, and are therefore potentially involved in explaining the modern rockfall activity at Hungerli Valley. The slope angle of the rockfall (low eigenvalue of −0.04) has low correlation to the other variables. We suggest this is due to the steep nature of the rockwalls, such that threshold slope conditions and sensitivity to triggers exist over much of the rockwall 68,69 . Consequently, at the scale of individual rockwalls, the external drivers (paraglacial stress release and debuttressing, and periglacial processes), all of which are elevation dependent, may be more important than slope form (i.e. slope oversteepening from glacial erosion) at explaining the local distribution of rockfall 5,70,71 . Rockfall volume (−0.20) relates to PC2. Rockfall volume is likely to be influenced by differences in the fractures spacings and patterns of the rock masses present in the valley, which while not explicitly included in the PCA, are captured by variation in elevation and possibly slope angle.  Spatial distribution of rockfall driving processes. Here we refine and analyse the spatial distribution of the three most correlated PCA variables; frost cracking, MARST (permafrost) and deglaciation age (paraglacial adjustment). Whereas for the PCA we treated elevation as a proxy for frost cracking, here we modelled frost cracking by applying a thermo-mechanical model 72 to measured rock surface temperatures 32 . Modelled frost cracking affects all rockwall locations in the Hungerli Valley, but is most intense within the three highest permafrost-affected rockwalls (RW1-3; Supplementary Fig. 5), and varies by rock type/strength within each elevation band ( Supplementary Fig. 7). Rock type ranges from weaker schisty quartz slate with aplite inclusions in the Rothorn cirque to stronger schisty quartz slate with inclusions of amphibolite and aplite in the Hungerlihorli (Supplementary Figs. 1 and 2) 73 . Regardless of the influence of lithological variations (Fig. 4b), a strong (non-monotonic) relationship was persistent between frost cracking and elevation. Our model results suggest frost cracking increases with elevation up until 3000 m, thereafter reducing at higher elevations (Fig. 4b). The modelled penetration depth of frost cracking is mostly shallower than 0.5 to 2 m, and exhibits increasing depth at higher elevations, up until 3000 m ( Supplementary Fig. 6). The modelled frost cracking pattern corresponds to the measured fracture spacing 32 , with more intense frost cracking generating more near-surface fracturing. As potential rockfall size is limited by the spacing of fractures, areas of more intense frost cracking will likely produce more numerous but smaller rockfalls, and be unlikely to trigger large events (>10 m 3 ).
Modelled MARST (a proxy for permafrost) is inversely related to elevation in the Hungerli Valley (Fig. 4d), such that colder MARST and permafrost likelihood increase at higher elevations. Rockwall with MARST below −3°C is located at a mean elevation of 3113 m (1 st Q. = 2989 m; 3 rd Q. = 3190 m), whereas rockwall with MARST between −1 to 0°C is located at a mean of 2901 m (1 st Q. = 2765 m; 3 rd Q. = 3034 m) (Fig. 2d). Permafrostfree areas (taken to have positive MARST) exist mostly below 2760 m, whereas above 2800 m, 75% of the rockwall area is frozen, and above 3100 m nearly the entire rockwall is affected by permafrost (Fig. 2). While MARST provides an indication of permafrost presence 66,74 , it cannot be used to assess the depth of seasonal warming and thaw, which influence rockfalls by decreasing shear resistance 21,[75][76][77] or increasing hydrostatic pressure 21 . The depth of warming and thaw depend additionally on seasonal temperature fluctuations and insulating snow cover 78 , which are challenging to model at a landscape scale. Consequently, we use MARST only as a proxy for permafrost presence, which we in turn use to identify rockwalls that are likely to have permafrost-related instability processes present (without modelling the processes themselves).
The deglaciation history since the LGM has resulted in a general trend for newly exposed rockwalls to occur at higher elevations over time (Fig. 4c). The general trend deviates because glacier retreat can involve thinning (vertical retreat) and terminus recession (planform retreat), both of which expose new rockwall but at different elevations (as evidenced by historic photos; Supplementary Fig. 3 and Fig. 2a). Between the LGM and Egesen, the thickness of the valley glacier decreased, and newly-exposed rockwalls (mean Short-term rockwall erosion rates. To quantify short-term rockwall erosion in the Hungerli Valley, we converted our contemporary rockfall data (volumes, m 3 ) into annual rockwall erosion rates (mm a −1 ) by dividing by the area of scanned rockwall. We then explored how erosion rate varies with the patterns of frost cracking, permafrost and paraglacial processes described above. The mean rockwall erosion rate of the entire scanned area was 0.86 mm a −1 between 2016 and 2019, but varies with elevation (assessed for 100 m elevation bands; Fig. 4f). Between 2600 and 3000 m the mean rockwall erosion rates increased with elevation, from 0.02 mm a −1 at the 2600-2700 m elevation band to 2.04 mm a −1 at the 2900-3000 m band. At higher elevations, the rates became smaller (1.52 mm a −1 at 3000-3100 m; 1.73 mm a −1 at 3100-3200 m, and 0.02 mm a −1 above 3200 m). The elevation bands with the greatest erosion rates, between 2900 and 3200 m (Fig. 4f), correspond to elevations with the most intense frost cracking (Fig. 4b). While frostcracking rates are known to vary with elevation 32,72,79 , hitherto frost-cracking rates had not been quantitatively and directly linked to rockwall erosion rates. Previous research suggested that frost cracking is positively related to catchment-wide erosion but inversely related to elevation 80 . This contrasts with our observed positive relationship between rockwall erosion rates and elevation. Catchment-wide erosion rates are measured from stream sediments that integrate multiple erosion processes throughout a catchment 81 . Moreover, rockwall erosion products (i.e. coarse rockfall debris) tend to be stored (e.g. in talus) with poor downstream connectivity. Consequently, and especially because rockwalls are underrepresented according to hypsometric curves 3,7 , the contribution of rockwall erosion to catchment-wide erosion and the relationships between rockwall erosion and its drivers are poorly captured by catchment-wide erosion rate measurement approaches.
Rockwall erosion rates in the Hungerli Valley also vary with deglaciation age, generally increasing from the oldest to the most recently deglaciated rockwalls. Ice-free rockwalls (on the Hungerlihorli) and rockwalls deglaciated between the LGM and Egesen had erosion rates of 0:25 þ0:14 À0:10 and 0:04 þ0:01 À0:01 mm a −1 , respectively. Allowing for the possibility of a cold-based glacier occupying the Hungerlihorli (i.e. reducing the area of exposed rockwall) between LGM and Egesen, the erosion rate increases to 0:13 þ0:06 À0:0 mm a −1 . Between the Egesen and 1882, a large part of the cirque became glacier-free and 30% of all recorded rockfalls, including the largest (159 m 3 ), occurred within this area. These rockfalls contributed to an erosion rate of 0:79 þ0:13 À0:20 mm a −1 . For rockwalls deglaciated between 1882 and 1929 there was an erosion rate of 0:24 þ0:13 À0:13 mm a −1 , which increases to a rate of 0:92 þ0:88 À0:88 mm a −1 for those deglaciated between 1929 and 1964. The greatest rate, of 5:69 þ0:44 À0:44 mm a −1 , was for rockwalls deglaciated between 1964 and 2019 (Fig. 4g). In general, rockfalls were prevalent in proximity to the extant glacier, with 77 rockfalls or 29% of all measured events, including four >10 m 3 sourced from rockwalls deglaciated between 1964 and 2019. Enhanced rockfall erosion in proximity to a retreating cirque glacier margin has been observed elsewhere 44 for which paraglacial processes such as thermo-mechanical stresses 13 , unloading stresses amplified by cirque glacier erosion 12,82 , or increased frost cracking within a bergschrund 83 may play a role. While we do not know which processes operate at the Hungerli Valley, the rockwall erosion rate patterns, at least partly, appear to reflect paraglacial adjustment 9,84 . Moreover, while some studies have suggested paraglacial rockwall adjustment (i.e. stress release and fracture propagation) operates at a millennial timescale 12,15 , our data show that adjustment can be more rapid (i.e. over decadal timescales). The older rockwalls in the Hungerli Valley may still be experiencing gradual paraglacial adjustment, but the rates of adjustment have slowed over time, as predicted by paraglacial exhaustion curves 9,84 .
Long-term rockwall erosion rates. To evaluate rockwall erosion rates within a longer temporal framework, we used geophysics to measure the volume of three talus accumulations in the Hungerli Valley. Talus slope-1 (TS1)-within the cirque and below a rockwall source area with elevation of 2900-3000 m (Fig. 2a)integrates 55 years of rockfall accumulation, giving a rate of 50:7 þ10:9 À10:3 mm a -1 . This erosion rate is one order of magnitude higher than those measured using the contemporary (i.e. 2016-2019) rockfall data from the same (2900-3000 m) elevation band (2.04 mm a −1 ; Figs. 4f and 5f) or from rockwalls of equal deglaciation age (5:69 þ0:44 À0:44 mm a -1 ; Fig. 4g). Talus slope-2 (TS2) and Talus slope-3 (TS3) are located at the base of Hungerlihorli peak, an area last glacierized during the Egesen stage (Younger Dryas), and integrate 10 to 12 ka of rockfall activity. The contributing rockwalls are located between 2600 and 2891 m for TS3 and between 2620 and 2891 m for TS2, and their talus accumulations suggest rockwall erosion rates of 1:0 þ0:6 À0:4 mm a −1 (TS3) to 1:2 þ0:6 À0:5 mm a −1 (TS2). The erosion rates are one to two orders of magnitude higher than those measured using the contemporary rockfall (i.e. laser scanning) data from the equivalent elevation bands (0.02 to 0.31 mm a −1 ; Figs. 4f and 5f). It is probable that the rates from long-term measurements exceed those of the short-term measurements due to integrating an aggressive period of paraglacial adjustment (expected soon after glacier retreat), and more intense permafrost thaw and frost cracking during colder phases (e.g. LIA) of the Holocene. The long-term measurements possibly also integrate episodic co-seismic rockfall activity, with more than 36 M w 4.0 and larger events recorded in Switzerland since 250 AD 86 and possibly increased seismic activity following deglaciation like found elsewhere 17,18 .
Elevation-dependent paraglacial and periglacial processes and their implications for rockwall erosion patterns in the European Alps. The short-term and long-term erosion rates in the Hungerli Valley both share the same positive relationship between elevation and erosion rates that was found for data compiled for the entire European Alps (Fig. 5a), especially at elevations above 2500 m (Fig. 5b). Our analyses from the Hungerli Valley suggest that this relationship is driven by frost cracking, permafrost occurrence and deglaciation age, all of which vary with elevation. Our modelled frost cracking magnitudes are elevation-dependent (Fig. 4b) and correspond to erosion rates increasing with elevation (Fig. 5c) up until a threshold elevation. Above this threshold, frost cracking decreases with increasing elevation 72,79 due to a decreasing availability of water (essential for frost cracking) at higher elevations 72,87 . This decreasing pattern is reflected in our frost cracking results, with reduced frost cracking at 3150 m (Fig. 4b), but the decrease will be more apparent for rockwalls that reach elevations higher than those of the Hungerli Valley. It is likely that frost cracking will be relatively minor on the highest summits of the Alps in the present climate, but will increase with further warming.
When we substitute elevation for MARST (i.e. permafrost), we see that permafrost presence contributes to greater erosion rates at higher elevations at our research area (Fig. 5e). Based on the European Alps permafrost model, most rockwalls above~2800 m are potentially affected by permafrost 38 , and therefore are likely to be experiencing elevated rockwall erosion rates, as presumed by previous high-elevation studies [41][42][43]45,46 . The development of regional scale models that can quantify the magnitude of activelayer thaw and permafrost warming effects (in addition to permafrost presence) will help to further evaluate the patterns we observe.
Substituting elevation for deglaciation age (Fig. 4c) reveals that short-term erosion rates increase with elevation in the Hungerli Valley (Fig. 5d), suggesting an elevation dependency to paraglacial adjustment. Over 80% of glaciers in the Alps are very small and located above 2500 m. These glaciers are very sensitive to climatic changes 37 and are retreating rapidly (in length but also laterally, especially within cirques), as observed for the Rothorn Glacier (Fig. 2). Consequently, our catchment scale study suggests that paraglacial processes will be responsible for a short-term increase (i.e. rapid adjustment) in rockwall erosion, especially at elevations above 2500 m.
Permafrost, frost cracking, and glacier retreat do not operate in isolation, but overlap and interact. Like for most alpine valleys, rockwalls deglaciated since the LIA in our research area (Fig. 4c) are located at elevations with permafrost (Fig. 4d). When warm-based glaciers retreat within permafrost terrain, rockwalls lose the insulating effect of the ice 88,89 , allowing permafrost to encroach into the newly exposed rockwall 89 , potentially increasing rockwall erosion 23 until the climate has sufficiently warmed to degrade all permafrost at those elevations. In a similar way, glacier retreat can introduce new thermo-mechanical stresses by exposing rockwalls to ambient temperature fluctuations 13 . Glacier retreat can also expose large bedrock fractures lengthened by glacier stresses 12,14,15 . As frost cracking is sensitive to crack length 32,72,87 and larger cracks freeze first 90 , the exposure of large cracks can facilitate more efficient frost cracking 32 . Permafrost presence can also amplify frost cracking, through bidirectional freezing 22 , which produces the highest rates in frost cracking models 32 . As a consequence of these interactions, climate warming has a profound, but non-uniform, effect on the breakdown of rockwalls by shifting the distribution and intensity of multiple paraglacial and periglacial processes. Climate warming in the European Alps is causing the retreat, and eventual disappearance, of very small glaciers 37 , which will drive paraglacial rockwall adjustment and the temporary growth-but long-term warming and degradation-of permafrost 33 . That warming and retreat will intensify frost cracking at higher elevations and in proximity to retreating glaciers, respectively, whereas at lower elevations the warming will reduce frost cracking 32 . Although lower elevations will experience a reduction of periglacial and paraglacial processes, other climate-dependent weathering processes 19 such as chemical and thermal weathering 91-93 may replace them as the main drivers of rockwall erosion 94 . Over time, the areas of most intense periglacial and paraglacial processes and resulting erosion will shift to higher elevations, in locations where higher topography exists.

Methods
Rockwall erosion rates in the European Alps. We compiled existing rockwall erosion rate data for the European Alps (Supplementary Table 1) and classified them into short-term (≤10 years), mid-term (>10 to ≤2500 years) and long-term (>2500 years) erosion rates. Short-term erosion rates were derived from rockfall collectors 41 or terrestrial laserscanning campaigns 40,42,43,[45][46][47] and integrate several years of rockfall activity (Fig. 1d). Due to the short integration time, these erosion rates are unlikely to include low-frequency high-magnitude events; however, if they do then the rates can be grossly over-estimated. Mid-term erosion rates (integrating 10 to 2500 years) were derived from lichenometric measurements 50 , averaged rates from talus-derived rock glacier deposits 48 or up-scaled rates using rockfall frequency-magnitude relationships 43,49 . Long-term rockwall erosion rates (integrating > 2500 years) were provided from geophysical imaging [51][52][53][54] or morphological estimation 40,58 of talus slope accumulations or rock glacier deposits 57 , and from cosmogenic nuclide exposure age dating 31,55,56 . Mid-to long-term rates potentially capture paraglacial and periglacial processes from pre-historical glacier recession and associated climate conditions (e.g. LIA or Younger Dryas; Fig. 1d). We compiled the mean, and the lower and upper range of the erosion rates to reflect the variability of erosion within catchments, and evaluated the relationship between erosion rate and elevation (Fig. 1b). If rockwall elevations were not reported in the original studies, they were derived from topographic maps or Google Earth.
Reconstruction of glacier retreat history of the Hungerli Valley. We used an orthophoto from 2017, historical photos, moraine locations and a regional study on LGM glacier thickness 95 to reconstruct the deglaciation history of the Hungerli Valley. To assess the locations and elevations of newly exposed rockwall as a result of glacier retreat since the LGM 96 , we used a regional study that suggests an ice surface elevation of 2800 m 95 for the main trunk glacier in proximity to Hungerli Valley at the LGM. Based on the ice surface isoline, the glaciated areas at the LGM were coherently enlarged to the cirque and firn basin headwalls 97 . For prehistorical rockwall exposure, we reconstructed the Egesen glaciation from mapped moraines ( Supplementary Fig. 3e), previously correlated to Egesen age (11 ± 1 ka 96 ; equivalent to Younger Dryas) 63,98 . We projected the moraine height to the rockwall (Fig. 2a) to approximate the ice surface at the rockwalls.
These approximate ice surface reconstructions suggest that the Rothorn cirque was fully glacierized during LGM and Egesen, while the upper Hungerlihorli peak (>2800 m) remained ice-free. While our reconstruction approach provides only a crude approximation, we consider it acceptable for the purpose intended, and we account for uncertainty by applying a ± 50 m vertical buffer to our Egesen and LGM glacier extents, which is propagated into the calculated erosion rates. Reconstructions based on trimlines ignore small, cold-based glaciers that might have been present on higher parts of the rockwalls, and to examine the impact of this we included a scenario with a cold-based glacier on the Rothorn during the LGM.
Historical photos were used to reconstruct the exposure of new rockwall from the 19 th century onwards. Historic photos were available for 1964 99 , 1929 and 1882 ( Supplementary Fig. 3a-c), and swissTopo provided a georeferenced orthophoto from summer 2017 (Fig. 2b). We reconstructed the glaciated area of the Rothorn cirque by using distinct landforms, ridges and mountain peaks, with an estimated vertical error of ±5%.
Rock temperature measurements. To model frost cracking and mean annual rock surface temperature (MARST), we measured rock surface temperature (RST) at five rockwalls along a topographic gradient from 2580 to 3148 m (Fig. 2) using Maxim iButton DS1922L temperature loggers installed in 0.1 m deep boreholes, similar to previous studies 32,78,91,100 . The loggers recorded RST at 3 h intervals ( Supplementary Fig. 4) between 1 September 2016 and 31 August 2019 at RW2 and RW5, and between 29 August 2017 and 28 August 2019 at RW1. MARST was calculated to compare temperature logger locations and for permafrost modelling (Supplementary Table 2). The rockwall thermal regime is affected by snow cover, with a significant influence on permafrost distribution 78,100-102 and frost cracking 32,103 . We examined the daily standard deviation in RST as a measure of snow cover, and applied a uniform threshold of <0.5°C for positive and negative RST following established practice 78,100 (Supplementary Note 1).
Rock temperature modelling. We propagated measured rock surface temperature into the rockwalls by conductive heat transfer with latent heat effects at the phase transition between water and ice 32,103 . Rock temperatures T over depth z and time t in one dimension were calculated by conservation of heat and Fourier´s Law 104 : where Q h is vertical heat flux, ρ r is rock density, c is specific heat capacity, and κ is thermal conductivity. Following the approach of previous studies 32,103,105 , we used a temperature-dependent volumetric specific heat capacity C and thermal conductivity κ to consider latent heat L effects for the phase changing envelope DT, which was set between −1 and 0°C: : where the subscripts f and u stand for frozen and unfrozen, ρ b is bulk (rock, water, ice) density, W reflects the total water content of the substrate as a fraction of mass, and W u is the remaining unfrozen water content below T f = −1°C. W equals the porosity for unfrozen rock. We set W to 3% following previous rock temperature models 32,78,89,106 and assumed a fully saturated rockwall, consistent with previous investigations of rock moisture measurements 107,108 . For unfrozen water content W u , we set a value 5% 32,103 . Thermal conductivity κ and specific heat capacity c were geometrically weighted in the frozen and unfrozen region of the rock mass: where i, w, r are the subscripts for ice, water, rock, and ν is the fractional content. The geology of the Hungerli Valley is predominantly paragneiss from the Mischabel nappe, consisting of schisty quartz slate with thin bands of aplite at the Rothorn cirque and intersections of amphibolite and aplite at the Hungerlihorli 109 ( Supplementary Figs. 1 and 2). For the calculations, we used previously measured rock density data (Supplementary Table 3) and took published thermal conductivity λ and heat capacity c 110 values for these rock types. To determine the temperature depth profile in each rockwall, we used our logged temperatures as rock surface temperatures, and applied Eqs. 1-4 with a spatial resolution Δz of 0.1 m. To evaluate the start conditions, we used a spin-up of 4 years on each temperature logger until steady-state rock temperature conditions were achieved 32,78,105 . Rock temperature modelling revealed permafrost occurrence with 2 to 3 m active layer thawing at RW1-3, while seasonal freezing penetrated 2 to 4 m at RW4-5 ( Supplementary Fig. 5).
Frost cracking modelling. To model frost cracking, we used our modelled rock temperature profiles to drive a frost cracking model (Supplementary Note 2). We applied a model that incorporates rock strength (Supplementary Table 3) and which assumes a porosity change to occur alongside frost cracking 72 . The model determines an upper temperature limit ΔT c for segregation ice growth depending on rock strength: where T m is the bulk melting temperature, ρ i ice density, L latent heat, P c the critical ice pressure, K C the critical fracture toughness, and x i the crack radius. We incorporated rock strength by using critical fracture toughness previously determined for each lithology 73 . The lower limit of ice segregation is controlled by water availability, which is determined by permeability: where k pf is the frozen and k p0 the unfrozen permeability, ΔT f is the pore freezing point, α is a power law exponent, and ΔT is the temperature below the bulk freezing point. Water migration in the frozen fringe was calculated based on the unfrozen permeability of the rock mass, with a value of 10 −18 m 2 adopted for all rockwalls (based on the permeability of schisty quartz slate 111 ). While permeability is likely heterogeneous within a single rockwall, and will vary slightly between the amphibolite and aplite lithologies, we consider it reasonable to assume a uniform value, consistent with existing model approaches 72,83 . Sensitivity tests on permeability, initial crack length and fracture toughness demonstrated that while these factors affect the magnitude of frost cracking, they do not change the pattern of affected rock mass depth 32 .
To analyse the dependence of frost cracking on elevation, we calculated the total amount of permanent rock expansion (λ) induced by frost cracking at each rockwall 72 . Total rock expansion is a sum of the frost cracking intensity at different depths within the rockwall and is calculated here as: The diffusivity D is where α is a power law exponent and μ is the water viscosity 72 . Diffusivity is used as a measure of the propensity for frost-induced porosity changes. We calculated the total expansion separately for each year and lithology, to calculate a minimum, mean and maximum annual rock expansion total over the monitoring period ( Supplementary Fig. 7). We consider the total rock expansion to represent the amount of (depth-integrated) frost cracking taking place at each rockwall.
Rockwall permafrost model. We used a multiple linear regression model fitted with multi-annual rock surface temperature (RST) measurements 67 to model rock permafrost distribution in the Hungerli Valley, taking into account the effects of snow cover limiting incoming solar radiation 32 . We simulated the potential incoming solar radiation (PISR) for each 2 m digital elevation model (DEM) raster cell using the area solar calculator of ESRI ArcGIS 10.7.1 for a snow-free period 74,112,113 from June to October in hourly steps assuming an atmosphere with no diffusivity 67,105,114,115 . The 'long-term' average length of the snow-free period was derived from a 20-year meteo station dataset 116 , located at Oberer Stelli Glacier in the Matter Valley. The station is 2-3 km SE of our study site at 2910 m elevation, and the snow duration at the meteo station corresponded well to snow cover in the Hungerli Valley for the 3-year period observed by our rock surface temperature measurements ( Supplementary Fig. 8). To develop a statistical relationship between mean annual rock surface temperature (MARST) and the mean annual air temperature (MAAT) and PISR, we used temperature loggers at north and south facing slope aspects in the Hungerli Valley (two south facing loggers, RW1-S and RW2-S, and five north-facing loggers, RW1-5; Supplementary Table 4). We calculated MARST for each raster cell as: with a, b and c representing calculated coefficients. For every logger location, we calculated the MAAT for the respective elevation using a temperature lapse rate of 6°C km −1 derived from a comparison of daily air temperature between meteo station Oberer Stelligletscher at 2910 m and a near-by meteo station, Grächen at 1605 m 117 .
We tested the quality of the relationship using r 2 , mean absolute error (MSE) and root mean square error (RMS) (Supplementary Table 5). Similar to previous approaches, we applied a 30-year MAAT for the period 1990 to 2019 (Supplementary Table 4) to Eq. (9) for the rockwall area to incorporate the slow response time of rock permafrost 67,114 . We interpret a MARST ≤ 0°C as permafrost 74 .
Laserscanning and rockfall detection. We scanned the rock faces in our study area from the opposing side of the valley floor, with a total of six positions (Fig. 2a), providing comprehensive coverage of all rockwalls. Due to the steep nature of the rockwalls, the oblique look-angle of our terrestrial laser scanning resulted in negligible scan shadowing (overcoming a limitation of aerial LiDAR scanning of steep rockwalls 118 ). Only minor, gently sloping ledges, of little importance for rockfall generation 119 , were missed (confirmed by aerial LiDAR inspection). Laserscanning took place over four consecutive summers from 2016 to 2019 using Riegl VZ400 or Riegl VZ2000i scanners with identical accuracy (Supplementary Table 6). To produce a reference point cloud, we registered the 2017 scan positions by manually selecting several pairs of pseudo-homologous points. The point clouds were then aligned using an iterative-closest-point (ICP) algorithm 120,121 implemented as the Multi Station Adjustment (MSA) in RiScan Pro 2.9. The basic assumption for this registration approach is that topographic change (i.e. rockfall) between scan epochs is negligibly small compared to the total scanned surface area used for point cloud alignment. To minimise registration errors, we excluded any non-rockwall surfaces (e.g. debris slopes). We applied the MSA algorithm to every subsequent scan from 2016, and 2018 to 2019 and registered each point cloud to our reference scan from 2017. The registration error lies between 2.7 and 5.0 cm, similar to that reported by others 43,46,47,122,123 . An octree filter implemented in RiScan Pro was used to adapt all point clouds to a spatially continuous resolution of 0.1 m. We assumed that a rockfall event produces a distance difference between two point clouds, which we identified using the Multiscale Model to Model Cloud Comparison (M3C2) algorithm 124 implemented in Cloud Compare 2.9.1, similar to previous studies 47,125,126 . The algorithm parameterizes surface roughness and alignment uncertainty to estimate a confidence interval, with a rockfall or movement recognised when the confidence interval is exceeded 127,128 . In our analyses, we used the significant change option for determining rockfall with a projection scale of 1 m and a maximum depth of 15 m. Every identified potential rockfall event was analysed manually and visually, with its volume estimated using the cutfill volume tool in RiScan Pro. Due to the annual resolution of our data collection (i.e. long time between scans), it is likely that many 'events' assumed to be single rockfalls could have been multiple smaller rockfalls 129 , which would affect our magnitude-frequency analyses. However, this issue has little bearing on the overall spatial patterns of our annual rockwall erosion rates, for which an annual scanning frequency is deemed suitable 126 .
Statistical analysis of influencing factors. Environmental variables are often auto-correlated and a principal component analysis (PCA) can be used to find an axis or dimension that summarizes this redundancy 130 . We applied a PCA using R software packages FactoMineR and ggplot to evaluate the relationships between rockfall variables similar to previous approaches 131,132 . For each rockfall event, we extracted (i) elevation as a proxy for frost cracking, (ii) MARST as a proxy for permafrost, (iii) deglaciation age as a proxy for paraglacial adjustment, (iv) slope angle as a proxy for shear stress, and (v) rockfall volume. In using deglaciation age as a proxy for paraglacial process intensity we assumed that paraglacial processes were most intense at the start of deglaciation and decreased thereafter in the form of an exhaustion curve 9,84 . We used the rockfall elevation as a proxy for frost cracking, which at the Hungerli Valley has been shown to strongly scale with elevation 32 despite variability in lithology and other local influences 73,87 . For each rockfall event, we extracted slope angle and rockfall volume from our laserscanning-derived DEM. All parameters were standardised using where x is the mean and σ is the standard deviation 130 . To standardise slope angle, we converted degree angles to radians. We kept principal components with values >1 (PC1) or close to 1 (PC2; Supplementary Table 7 and Fig. 3b).
Derivation of short-term rockwall erosion rates. For evaluating the influence of elevation on rockwall erosion rates, we derived rockwall area within 100 m elevation bands (Supplementary Table 8). Based on our reconstructed glacier distribution, we calculated rockwall erosion rates for every deglaciation step (i.e. successively exposed rock wall area) incorporating a glacier extent uncertainty of ±5% for 2017, 1964, 1929 and 1882 glacier extents and an uncertainty of ±50 m for Egesen and LGM glacier extents (Supplementary Table 9). We classified MARST into <−3°C, −3 to −2°C, −2 to −1°C, −1 to 0°C and >0°C classes, similar to previous permafrost model approaches 74 , and calculated the area and rockwall erosion rate for every MARST class (Supplementary Table 10). Area was calculated in Cloud Compare using the Poisson Surface Reconstruction tool 133,134 to produce a mesh (with a relatively high octree depth of 10 to 12). The mesh was adjusted using the output density to optimise it for the input cloud, to capture the full complexity of the surface area, without oversmoothing data gaps and patchy areas of the point cloud.
Mid-and long-term rockwall erosion calculation using talus slopes. We selected only talus slopes that have not been reworked by rock glaciers for our analysis in the Hungerli Valley. We applied electric resistivity tomography (ERT) and seismic refraction tomography (SRT) in longitudinal transects from the rockwall to the toe of three talus slopes (Supplementary Note 3) to quantify the sediment thickness. For topographic corrections, we extracted elevations along each transect from the swissTopo Alti3D digital elevation model. ERT data were acquired using a Wenner-Schlummberger approach and were processed using Res2DInv with a resolution of half of the electrode spacing (Supplementary Table 11). The root mean square error (RMS) was used to provide an estimation of model quality, and the values were within the ranges found in previous talus slope studies 51,[135][136][137][138] . Seismic refraction signals were triggered by sledgehammer shots, recorded using a Geometrics Geode with 24 Geophones, and processed using ReflexW 9.1.3 (Supplementary Table 11). First-arrivals were picked manually 139,140 to derive the travel-time of the refracted seismic waves, which were analysed to identify seismic layers [139][140][141] . Layer properties were used to adjust initial starting models and gradients. Sediment thickness was interpreted from one-dimensional profiles of electrical resistivity and seismic velocity every 12.5 m ( Supplementary  Fig. 9). Talus slopes 2 and 3 indicated the possible presence of a basal moraine, which would result in smaller talus volumes. We accounted for this uncertainty by using two different basal talus surfaces in the volume calculations (with and without moraine). Talus slope volume (V t ) was calculated by subtracting the lower subsurface boundary from the talus ground surface using Golden Software Surfer 17. The rockwall erosion rate E r was calculated as: where ρ r is the density of the rockwall (2800 kg m −3 , Supplementary Table 3), ρ t is talus density, A r is the rockwall surface area providing rockfall material, and t is the time period of talus accumulation. Talus density varies vertically from 2000 kg m −3 near the surface to 2500 kg m −3 at the base 142 . We incorporated the entire range of talus density into our erosion rate calculation to avoid bias from densification 61 if using only a single value. Talus age was derived from glacier reconstruction and ranged from 55 a (Talus slope-1) to 11 ± 1 ka (Talus slope-2 and -3). Rockwall surface area and talus slope surface area were both determined from the DEM with an assumed mapping uncertainty of ±5%, which we included, along with potential moraine presence, in our erosion rate calculation (Supplementary Table 12). Mid-to long-term erosion rates integrate changes in seismic activity 86 or climatic phases 60 and the influence should be similar in a single mountain range as the European Alps.
Code availability