Last-glacial-cycle glacier erosion potential in the Alps

The glacial landscape of the Alps has fascinated generations of explorers, artists, mountaineers, and scientists with its diversity, including erosional features of all scales from high-mountain cirques to steep glacial valleys and large overdeepened basins. Using previous glacier modelling results and empirical inferences of bedrock erosion under modern glaciers, we compute a distribution of potential glacier erosion in the Alps over the last glacial cycle from 120 000 years ago to the present. Despite large uncertainties pertaining to the climate history of the Alps and unconstrained glacier erosion processes, the resulting modelled patterns of glacier erosion include persistent features. The cumulative imprint of the last glacial cycle shows a very strong localization of erosion potential with local maxima at the mouths of major Alpine valleys and some other upstream sections where glaciers are modelled to have flowed with the highest velocity. The potential erosion rates vary significantly through the glacial cycle but show paradoxically little relation to the total glacier volume. Phases of glacier advance and maximum extension see a localization of rapid potential erosion rates at low elevation, while glacier erosion at higher elevation is modelled to date from phases of less extensive glaciation. The modelled erosion rates peak during deglaciation phases, when frontal retreat results in steeper glacier surface slopes, implying that climatic conditions that result in rapid glacier erosion might be quite transient and specific. Our results depict the Alpine glacier erosion landscape as a time-transgressive patchwork, with different parts of the range corresponding to different glaciation stages and time periods.

1 Introduction

The glacial erosion landscape of the Alps has fascinated generations of explorers, artists, mountaineers, and scientists for centuries. Its cultural impact is indeed so far-reaching, that in English, a non-Alpine language, the adjective “alpine” with non-capital “a” is now casually used to describe an Alpine-like, glacially modified mountain landscape outside the Alps, while the proper noun “Alps” has been applied to nickname Alpine-like, glacier-eroded mountain ranges in Norway (Lyngsalpene), New Zealand (Southern Alps), Japan (, Nihon Arupusu), and elsewhere. Some mountain ranges are predominantly characterized by cirque glaciation (e.g. Uinta Mountains), glacial valleys (e.g. Putorana Plateau), or large-scale overdeepenings (e.g. Patagonia). But other regions, including the Alps (e.g. Penck, 1905), present a higher variety of glacier erosional landforms, whose implications for glacial history are yet to be understood. In parallel to such landscape diversity, the cosmogenic nuclide memory of bedrock erosion provides a more quantitative but equally varied picture of glacier erosion effectiveness both within and between glaciated regions (Jansen et al., 2019; Steinemann et al., 2020, 2021).

Glacially eroded topography has sometimes been used as a proxy for mapping palaeoglacier extent (e.g. Margold et al., 2011; Fu et al., 2012). Cold-based glaciers, however, have been observed to preserve landforms as fragile as sand beaches (Kleman, 1994), leading to the reinterpretation of complex glacial landscapes within a binary conceptual framework of cold and temperate basal thermal regimes (Kleman et al., 2008, 2010; Fabel et al., 2012; Fu et al., 2013). Nevertheless, it remains elusive whether glacial topography is a proxy for temperate ice-cover duration and why some glaciated regions show preserved periglacial blockfields topped by erratic boulders and others, including the Alps, do not (Wirsig et al., 2016; Seguinot et al., 2018).

On a larger scale, glacier erosion also likely governs, at least in part, the height of mountains (Egholm et al., 2009; Thomson et al., 2010) through the “glacier buzzsaw”. In turn, evidence suggests that increased mountain erosion rates coincided with global cooling and Pleistocene glaciations from roughly 2.6 million years before present (Herman and Champagnac, 2016). However, oceanic isotopic proxies for global weathering rates remained constant over this time period (Willenbring and Von Blanckenburg, 2010). Examination of contemporary erosion rates across different climates suggests that increased temperatures within glaciated regions lead to higher erosion rates (Koppes and Montgomery, 2009; Koppes et al., 2015; Fernandez et al., 2016). However, understanding the link between climate and glacier erosion in the sedimentary record remains difficult due to timescale biases (Ganti et al., 2016) and the non-linear relationship between atmospheric temperature and erosion (e.g. Anderson et al., 2012; Mariotti et al., 2021). Numerical glacier modelling may provide insights into such complexity.

Observing and quantifying long-term glacier erosion and sedimentation processes has been a challenge, and two general methods have been adopted to quantify erosion (Alley et al., 2019). Physically based models describe the quarrying process and abrasion of bedrock by debris-laden ice (e.g. Alley et al., 1997; Iverson, 2012; Beaud et al., 2014). Describing the processes physically provides a basis for understanding erosion processes (Hallet, 1979; Ugelvig et al., 2018) and yields insight into the formation of glacial landforms, such as tunnel valleys and eskers (Beaud et al., 2018; Hewitt and Creyts, 2019). Yet, these models prove difficult to implement in many cases due to the large number of poorly constrained parameters and processes, i.e. water pressure fluctuations and subglacial debris concentration (e.g. Hallet, 1979; Seguinot, 2008; Ugelvig et al., 2018).

Instead, empirical relationships between glacier sliding and erosion can represent erosional quantities well and can recreate important glacial landforms (e.g. Harbor et al., 1988; MacGregor et al., 2000). For instance, Humphrey and Raymond (1994) correlated temporal variations in ice velocity and suspended sediment load during a surge of Variegated Glacier to establish a linear erosion law. Koppes et al. (2015) quantified sediment yields in 15 Patagonian and Antarctic Peninsula fjords, concluding there was a predominant control of surface air temperature on rapid glacier erosion and calibrating a near-square law from glacier velocity to erosion rate. Herman et al. (2015) collected suspended sediment samples for 5 months in the outlet stream of Franz Josef Glacier, mapped their chemical composition to geological zones of distinctive glacier speed, and also concluded there was a near-square relationship between basal sliding and glacier erosion. Cook et al. (2020) assembled a global compilation of erosion rates for 38 glaciers, showing a predominant role of glacier sliding velocity over climate variables yet concluding there was a sub-linear relationship to erosion.

Here, the non-linear erosion law by Koppes et al. (2015) is applied to previously published model results (Seguinot et al., 2018) and the patterns of modelled erosion potential and cumulative last-glacial-cycle erosion potential are analysed. We examine the glaciological conditions that cause erosion and discuss the implications of these conditions on understanding the mark of ancient glaciers on Alpine topography and the relationship between climate and glacier erosion. Despite aggregated uncertainties in palaeoclimate, glacier flow, and glacier erosion processes, our results provide insights into the diversity of the Alpine glacial erosion landscape.

2 Methods

2.1 Ice sheet modelling

The ice sheet model set-up was presented in an earlier publication (Seguinot et al., 2018) and is briefly summarized here. The simulations were conducted using the Parallel Ice Sheet Model (PISM, development version e9d2d1f), an open-source, finite-difference, shallow-ice, and shallow-shelf model (the PISM authors, 2017). Our set-up includes temperature-dependent and water-content-dependent creep, pseudo-plastic basal sliding, bedrock isostatic deformation under the ice load, and a positive degree-day (PDD) surface mass balance model. The model is initialized with assumed present-day ice thickness and 1 km resolution bedrock topography from SRTM (Jarvis et al., 2008) and ice and bedrock temperature at 120 ka (1 ka is 1000 years before the present) and is run to the present. The ice sheet model physical parameters are listed in the earlier open-access publication (Seguinot et al., 2018), and the complete list of PISM parameters is stored in long-term, archived model output metadata for the present study (Seguinot, 2021a) and for the earlier publication (Seguinot, 2019, 2020).

Climate forcing is provided by a monthly climatology from interpolated observational data (WorldClim; Hijmans et al., 2005) and the European Centre for Medium-Range Weather Forecasts Interim Reanalysis (ERA-Interim; Dee et al., 2011), amended with temperature changes from the European Project for Ice Coring in Antarctica (EPICA; Jouzel et al., 2007) and time-dependent palaeo-precipitation reductions (Huybrechts, 2002). Coarser, 2 km resolution runs based on alternative palaeoclimate scenarios (Seguinot et al., 2018), forced by temperature changes from the Greenland Ice Core Project (GRIP; Dansgaard et al., 1993) and an oceanic sediment core from the Iberian margin (MD01-2444; Martrat et al., 2007), are analysed in the Discussion section (Sect. 4.1).

2.2 Erosion law

Modelled erosion rates, $\dot{e}$, are calculated from the modelled basal velocities, u_b, using an empirical erosion power law,

$$\begin{array}{}\text{(1)}& \dot{e}=K_{\mathrm{g}}{u}_{\mathrm{b}}^l,\end{array}$$

where $K_{\mathrm{g}}=\mathrm{5.2}\times {\mathrm{10}}^{-\mathrm{11}}$ m^1−l a^l−1 and l=2.34 are empirical constants calibrated on quantified fjord sediment yields and velocities from 15 Patagonian and Antarctic Peninsula tidewater glaciers (using the full dataset from Koppes et al., 2015, including outliers). We discuss the role of alternative values of K_g and l in the Discussion section (Sect. 4.2).

Our approach does not account for feedbacks of glacier erosion onto bedrock topography and then onto ice dynamics (e.g. Anderson et al., 2012). Furthermore, we do not account for the downstream advection of debris (e.g. Anderson and Anderson, 2016) or fluvial transport of subglacial sediment (Delaney et al., 2019). Instead, we assume that eroded material is instantly transported out of the system, thus neglecting its role in shielding the bedrock from glacier erosion in zones of temporary storage (Preusser et al., 2010). Nor do we account for differences in erosion effectiveness on different lithologies or erosion from subglacial and interglacial hydrological processes. For the aforementioned reasons, we refer to the above-computed rates, $\dot{e}$, as “potential erosion rates”.

Secondly, we compute the time-integrated glacier erosion rate,

$$\begin{array}{}\text{(2)}& e=K_{\mathrm{g}}\int u_{\mathrm{b}}^l\mathrm{d}t,\end{array}$$

and refer to it as the “cumulative erosion potential”. The integration is numerically approximated using a time step of 10 a in the main simulation and 100 a in the alternative climate scenarios runs. Note that in the above formula with an exponent, l, higher than 1, the cumulative erosion potential has a greater dependency on the basal ice dynamics, u_b, than the duration of glacier cover.

3 Results

3.1 Cumulative erosion

The modelled cumulative (time-integrated) glacial erosion potential (Fig. 1a) varies by several orders of magnitude from insignificant to 100 m scale erosion potential. Its spatial patterns show a very strong localization along the Alpine valleys, with local maxima occurring both at the Alpine gates, where ice flow transited from valley to piedmont glaciers, and further upstream, where the valley slopes increase. While mountain cirque glaciers and relevant erosion processes may not be captured by the model physics and horizontal resolution, high cumulative erosion potential also occurs near the valley heads. There is a general tendency for higher cumulative erosion in the north-western Alps where the input winter precipitation is higher (WorldClim; Fig. 1h in Seguinot et al., 2018) and the glacial relief more pronounced in the topography.

Figure 1 (a) Modelled cumulative (time-integrated) glacial erosion potential over the last glacial cycle and geomorphological reconstruction of Last Glacial Maximum Alpine glacier extent for comparison (solid blue line; Ehlers et al., 2011). The background map consists of the initial basal topography from SRTM (Jarvis et al., 2008) and Natural Earth data (Patterson and Kelso, 2017). (b) Modelled total ice volume in centimetres of sea-level equivalent (cm s.l.e., black), annual (domain-integrated) potential erosion volume (light brown), and its 1 ka running mean (dark brown). Shaded grey areas indicate the timing for MIS 2 (marine oxygen isotope stage) and 4 (Lisiecki and Raymo, 2005). Hatching marks periods with ice volume below 3 cm s.l.e. where glacier sliding may be affected by stress-balance approximations and the model horizontal resolution (see Sect. 4.4).

3.2 Temporal evolution

The modelled annual (domain-integrated) potential glacial erosion volume does not systematically correlate with the modelled total ice volume (Fig. 1b). Except for the onset and termination of the glacial cycle were ice volume is low and sliding processes are unlikely to be captured by the model resolution, a comparison between the modelled ice volume and the modelled annual erosion volume shows no single relation between these two quantities (Fig. 2). After significant ice volume has been reached, there is a general tendency for slower erosion during periods of extensive glaciation (Fig. 2). However, the domain-integrated erosion volume is modelled to be consistently higher (by a factor of 2 to 10) during periods of decreasing ice volume than during periods of increasing ice volume (Fig. 2).

Figure 2 Modelled annual (domain-integrated) potential erosion volume (light curves) and its 1 ka running mean (dark curves) in relation to the modelled total ice volume in centimetres of sea-level equivalent. Blue dashed lines indicate increasing ice volume and solid brown lines decreasing ice volume. Hatching indicates low confidence due to limited glacier volume.

3.3 Spatial migration

A closer look at Rhine Glacier, one of the palaeo ice sheet's largest outlets, reveals a spatial migration of the modelled rapid erosion areas. During stages of glacier advance and maximum extension, erosion is modelled to be under 1 mm a⁻¹ and restricted to the lower parts of the catchment, while much of the intra-montane Alps (modelled to be largely cold-based; Fig. 6c in Seguinot et al., 2018) experience insignificant erosion potential (Fig. 3a and b). The potential erosion rates both increase and propagate inwards during periods of glacier retreat (Fig. 3c and e).

Figure 3 (a–d) Modelled instantaneous potential erosion rate of Rhine Glacier for selected glacier advance and retreat ages and the final model state for topographic reference. (e) Interpolated instantaneous potential erosion rate along a Rhine Glacier transect (upper panels, dashed line) for the entire last glacial cycle.

The results found in the Rhine Glacier catchment can be generalized to the entire model domain by using (present-day) bedrock altitude as a proxy for along-flow distance (Fig. 4). Thus more generally, periods of modelled increasing and maximum ice volume correspond to lower values for elevation-aggregated potential erosion rates, with significant erosion potential restricted to lower elevations. In contrast, periods of modelled decreased ice volume correspond to higher local modelled erosion rates and more extensive rapid erosion potential in higher-elevation areas (Fig. 4).

Figure 4 (a) Potential erosion rate “hypsogram”, showing the geometric mean of modelled rates in 100 m elevation bands across the entire model domain and its time evolution. Hatching indicates elevation bands with fewer than a hundred grid cells (100 km²). (b) Distribution of model domain bedrock topography (grey bars) and glaciated topography (blue bars) in 100 m elevation bands and the corresponding cumulative potential erosion volume (brown line). (c) Same as Fig. 1b.

4 Discussion

4.1 Climate sensitivity

Different palaeoclimate forcing scenarios lead to different modelled glaciation histories (Seguinot et al., 2018) and thus different results in terms of glacial-cycle cumulative erosion potential (Fig. 5). Scenarios including higher (constant, present-day) input precipitation rates yield generally higher modelled ice discharge and thus higher erosion potential (Fig. 5a–c), while scenarios including palaeo-precipitation reductions (Fig. 5d–f), including our default high-resolution run (Fig. 1), yield a smaller erosion potential.

Figure 5 Modelled cumulative glacial erosion potential over the last glacial cycle using three different palaeo-temperature histories from (a, d) the Greenland Ice Core Project (GRIP; Dansgaard et al., 1993), (b, e) the European Project for Ice Coring in Antarctica (our default, EPICA; Jouzel et al., 2007), and (c, f) an oceanic sediment core from the Iberian margin (MD01-2444; Martrat et al., 2007), both (a–c) without and (d–f) with palaeo-precipitation reductions (Seguinot et al., 2018).

Figure 6 Modelled cumulative (time-integrated) glacial erosion potential for three different erosion laws published by (a) Koppes et al. (2015), $\mathrm{5.2}\times {\mathrm{10}}^{-\mathrm{8}}{u}_{\mathrm{b}}^{\mathrm{2.34}}$ (same as Fig. 1a but with a different colour scale); (b) Herman et al. (2015), $\mathrm{2.7}\times {\mathrm{10}}^{-\mathrm{7}}{u}_{\mathrm{b}}^{\mathrm{2.02}}$; (c) Humphrey and Raymond (1994), $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}{u}_{\mathrm{b}}^{\mathrm{1}}$; and (d) Cook et al. (2020), $\mathrm{1.665}\times {\mathrm{10}}^{-\mathrm{1}}{u}_{\mathrm{b}}^{\mathrm{0.6459}}$. (e–h) Corresponding modelled cumulative glacial erosion potential along a Rhine Glacier transect (upper panels, dashed line) and (i–l) corresponding modelled evolution of annual (domain-integrated) erosion volume (light brown) and 1 ka running mean (dark brown) in relation to modelled total ice volume in centimetres of sea-level equivalent (as in Fig. 2a).

The spatial pattern of modelled time-integrated erosion is nevertheless generally consistent across the palaeoclimate scenarios tested, despite variations in the glaciation maximum extent (Fig. 3 in Seguinot et al., 2018; Fig. 5). It should be noted, however, that all runs presented here show a systematic bias with excessive glacier cover in the Eastern Alps and a diminished glacier extent in the Western Alps (Fig. 1a; further discussed in Seguinot et al., 2018). Thus the modelled patterns of erosion potential certainly include a similar bias. However, uncertainties due to palaeoclimate appear minor in regard to the much larger uncertainties affecting the choice of erosion law.

4.2 Choice of erosion law

The choice of erosion law significantly impacts the results (Fig. 6). Our default, based on quantified sediment yields from Patagonian and Antarctic Peninsula tidewater glaciers ($\dot{e}=\mathrm{5.2}\times {\mathrm{10}}^{-\mathrm{8}}{u}_{\mathrm{b}}^{\mathrm{2.34}}$; Koppes et al., 2015) yields a moderate and strongly localized cumulative erosion potential (Figs. 1 and 6a). The erosion law based on measured suspended sediment load from Franz Joseph Glacier ($\dot{e}=\mathrm{2.7}\times {\mathrm{10}}^{-\mathrm{7}}{u}_{\mathrm{b}}^{\mathrm{2.02}}$; Herman et al., 2015) also yields a strongly localized, and yet much higher, modelled erosion potential (Fig. 6b). The erosion law based on measured suspended sediment load during a surge of Variegated Glacier ($\dot{e}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}{u}_{\mathrm{b}}^{\mathrm{1}}$, Humphrey and Raymond, 1994) results in similarly high values but a less localized pattern of cumulative erosion potential (Fig. 6c). Finally, the erosion law based on a global compilation of glacier velocities and erosion rates ($\dot{e}=\mathrm{1.665}\times {\mathrm{10}}^{-\mathrm{1}}{u}_{\mathrm{b}}^{\mathrm{0.6459}}$; Cook et al., 2020) yields an even flatter pattern of cumulative erosion potential (Fig. 6d).

With a total Pleistocene glacial relief on the order of 1 km (Preusser et al., 2011; Valla et al., 2011), a cumulative glacial erosion for the last glacial cycle on the order of 10 to 100 m can be expected. However, none of the tested erosion power laws fall within this range. Instead, the erosion law calibrated on tidewater glaciers (Koppes et al., 2015) yields cumulative erosion in the Rhine Valley on the order of metres, while the three erosion laws based on terrestrial glaciers (Humphrey and Raymond, 1994; Herman et al., 2015; Cook et al., 2020) result in kilometre-scale integrated erosion potential. During the Last Glacial Maximum and much of the last glacial cycle, Alpine palaeoglaciers were closer in size, slope (an important parameter as we argue in the next section), and climatic context to the present-day glaciers of Patagonia and the Antarctic Peninsula (Koppes et al., 2015) than to Franz Joseph Glacier (Herman et al., 2015) and many of the glaciers included in the global compilation by Cook et al. (2020). This may help to explain why the reality appears to fall in between the tested erosion laws.

That said, the magnitude of the modelled erosion rates can hardly be used as a criterion to rule out this or that erosion law for the glacial Alps. While sliding is modelled here by a pseudo-plastic law based on till dilatation under high water pressure (Tulaczyk et al., 2000), empirical erosion rules are thought to be a proxy for glacier abrasion (e.g. Hallet, 1979; and possibly plucking) of a glacier flowing over a hard bed. Current erosion rates in the Alps strongly vary depending on the substratum (Steinemann et al., 2020, 2021), and it is arguable that a layer of till and proglacial sediments may have armoured the glacier bed against erosion for at least parts of the glacial cycle. Thus, given the lack of data to constrain the erosion rule, it can not be excluded that a linear or sub-linear erosion law (e.g. Cook et al., 2020) would lead to a less localized pattern of erosion (Fig. 6c and d) and a stronger correlation between modelled ice volume and erosion (Fig. 6k and l). On the other extreme, an even more non-linear erosion law, not tested here, derived from tidewater glaciers but excluding two outliers ($\dot{e}=\mathrm{5.3}\times {\mathrm{10}}^{-\mathrm{9}}{u}_{\mathrm{b}}^{\mathrm{2.62}}$; Koppes et al., 2015) would result in an even more localized pattern of erosion potential.

The modelled pattern of erosion potential depends on PISM glacier physics and sliding model parameters. In particular, the pseudo-plastic sliding law exponent (q=0.25 in Seguinot et al., 2018) controls the sharpness of the transition between adherent and decoupled basal conditions. Recent ensemble validation of Antarctic Ice Sheet glacial-cycle simulations against geological and present-day observations (Albrecht et al., 2020a, b) supports a higher value of q=0.75 and thus a sliding law closer to linear, which would perhaps result in a smoother distribution of sliding velocity and erosion potential. Lateral stress gradients missing from the shallow-shelf approximation stress balance could also contribute to moderate sliding velocity in narrow troughs (Herman et al., 2011; Egholm et al., 2012a, b; Pedersen et al., 2014).

Because the intrinsically large variations in modelled glacier sliding velocities yield cumulative glacial-cycle erosion 2 to 3 orders of magnitude higher in the Alpine valleys than near the mountain tops (Fig. 6d), we argue that the modelled highly localized pattern of erosion potential is qualitatively robust. This pattern, as well as the lithologic variations in bedrock erodibility, may explain order-of-magnitude variation in glacier erosion rates observed in the cosmogenic nuclide record (Jansen et al., 2019; Steinemann et al., 2020, 2021).

4.3 Climate control on erosion

The modelled, domain-integrated glacier erosion volume does not increase together with ice volume. In contrast, periods of ice advance and maximum glaciation correspond to comparatively slow erosion potential (Figs. 1 and 2). This result may seem paradoxical and is modulated to a limited extent by the exponent, l, in the erosion rule (Fig. 6). Nevertheless, it may explain why some areas covered during glacial maxima only, such as some valleys on the southern side of the Alps, appear to have experienced only little glacial modification of their topography. During deglaciation periods, increased surface melt at low elevations results in a steepened topographic profile (Fig. 7a), as observed and further expected on currently retreating glaciers (Huss et al., 2010; Zekollari and Huybrechts, 2015). The buttress formed by the glacier foot against gravitational forces is reduced, causing an increase in the magnitude of the basal shear stress (basal drag) relative to the ice overburden pressure and an increase in the area where sliding and, thus, erosion can occur (Fig. 7b).

Figure 7 (a) Modelled basal and surface topographies interpolated along Rhine Glacier during glacier advance (36 ka, dashed blue line), maximum extent (24 ka, solid grey line), and retreat (16 ka, solid brown line, using the same transect as in Fig. 7). (b) Interpolated yield stress in the pseudo-plastic sliding law (dotted lines) and magnitude of the basal shear stress (basal drag, solid lines), both normalized against the ice overburden pressure. Basal drag values greater than yield stress result in glacier sliding (see Seguinot, 2014; Seguinot et al., 2016, for model physics). Some of the observed basal drag fluctuations in the upper part of the transect result from steep slopes on either side of the narrow valley.

This result is corroborated by a recent denudation record from the Mediterranean Alps, which includes an increase in glacier erosion in the Var Valley, roughly at 25 ka following an extended period of low denudation rates (Mariotti et al., 2021). While that increase in erosion pre-dates the increase observed in our study (ca. 15 ka; Fig. 2), it demonstrates that climatic conditions over which glacier erosion occurs are relatively specific and likely transient. Other field-based studies have evidenced increased glacier erosion and proglacial sedimentation during the current deglaciation period (e.g. Koppes and Montgomery, 2009; Micheletti and Lane, 2016; Lane et al., 2017; Bendixen et al., 2017). This increase has often been attributed to higher meltwater availability; promoting sediment transport (Delaney and Adhikari, 2020); enhanced glacier sliding (Herman et al., 2011); and hydrofracturing, i.e. plucking (Hallet, 1996; Ugelvig et al., 2018; Hildes et al., 2004).

But importantly, these hydrological processes are not included in our model, so deglaciation surface meltwater is assumed to instantly exit the glacier with no effect on glacier dynamics (cf. Werder et al., 2013; Iverson, 2012; Ugelvig et al., 2018). The enhanced glacier velocity during periods of ice retreat is here modelled to result from changes in glacier geometry alone (Fig. 7b). Coincidentally, the link between climate and erosion is also the topic of a recent, independent study using PISM (Lai and Anders, 2021). Observed feedbacks between surface meltwater and basal sliding, not included in the model, presumably further accelerated erosion during deglaciation (Herman et al., 2011). Furthermore, we omit the fluvial transport of sediment, which could well increase sediment discharge and erosion with increased melt, so long as additional sediment remains available (Delaney and Adhikari, 2020). Rapid glacier erosion can be expected to accompany the ongoing 21st-century deglaciation, especially for the ice sheets of Greenland and Antarctica where polythermal and cold-based glaciers are warming and sliding may increase (e.g. Moon et al., 2012; Mouginot et al., 2014; Overeem et al., 2017).

4.4 Age of the glacial landscape

In many places throughout the Alps, significant modelled cumulative erosion occurs far inside the Alpine valleys (Fig. 1). However, the potential erosion patterns experience spatial shifts through the glacial cycle (Fig. 3), so high-altitude rapid erosion is not modelled to date from expansive glaciation stages but instead from stages of intermediate glacier cover (Fig. 4; Barr et al., 2019).

Figure 8 Alpine glacial erosion landscape diversity. (a) Piedmont overdeepening of Lake Constance, ca. 10×50 km. (b) Glacial trough at Bout du Monde, ca. 1×10 km. (c) The unusually deep mountain cirque revealed by the current demise of Chüeboden Glacier, ca. 1×1 km.

For instance, low-elevation piedmont overdeepened basins such as Lake Constance (Fig. 8a), Lake Geneva, and Lake Maggiore are here modelled to be the product of extensive Last Glacial Maximum-like glaciations. Other modelling work suggests that features such as these likely form when water can access the glacier bed and encourage glacier sliding and, thus, erosion (Herman et al., 2011). The model herein does not include subglacial hydrology, yet the role of climate warming in increasing glacier erosion at the lower elevations remains evident (Fig. 4). Alpine terrain with elevations above 1000 m is modelled to experience virtually no erosion from 35 to 18 ka BP (Fig. 4), so many of the most spectacular, intermediate-elevation glacial valleys of the Alps such as Lauterbrunnental, Haslital (upper Aare), or Bout du Monde (upper Giffre, Fig. 8b) would date from periods of intermediate glacier extension.

The validity of the model results at high elevation is discussable. Crystalline massifs such as the Écrins, Gran Paradiso, Monte Rosa, Aare, Ötztal, and High Tauern massifs locally exhibit a strikingly high erosion potential. However, the computation of glacier flow velocities on such steep surfaces is strongly limited by the model horizontal resolution of 1 km, the shallow-ice glacier flow physics (Imhof et al., 2019), and PISM's current mass-conservation heuristics (Imhof, 2021, Appendix). Besides, bergschrund (rimaye) processes likely to dominate interglacial cirque erosion at such altitudes (Sanders et al., 2012) are not captured by the velocity-based glacier erosion power laws. Nevertheless, the lack of rapid erosion at high elevations for much of the modelled glacial cycle (Fig. 4) implies that the highest mountain cirques (e.g. Fig. 8c) would have formed over interglacial periods, when glaciers are confined to high terrain (Barr et al., 2017, 2019). Similar processes also occur on Himalayan glaciers and may limit erosion for the high-elevation portion of that mountain range during cold periods as well, failing to offset uplift rates (Harper and Humphrey, 2003). Over longer timescales, though, the vertical distribution of erosion rates also depends on the erosional modification of topography (Sternai et al., 2013).

Finally, the time-transgressive nature of the modelled glacial erosion rates hints at a new explanation for the Alpine glacier trim line and lack of cold-based glaciation evidence. Due to their geographic location in the mid-latitudes and their steep topographic gradient, the Alps appear to be a mountain range that hosted glaciers of various sizes throughout the varying climate of the Quaternary, resulting in a transgressive localization of glacier erosion throughout their elevational range, from the piedmont during glacial maxima to the highest cirques during interglacial periods. The Alpine trim line, in this case, would correspond neither to the upper reach of the Last Glacial Maximum glaciers (e.g. Kelly et al., 2004) nor to an englacial temperate-ice boundary (Coutterand, 2010; Seguinot et al., 2018) but to a time-transgressive upper limit of erosion from advancing and retreating glaciers on steep terrain.

5 Conclusions

Due to compounded uncertainties regarding palaeoclimate, glacier sliding, and erosion processes, our quantitative results are very likely inaccurate. However, we draw the following qualitative conclusions:

• The non-linear physics of glacier deformation and sliding combined with non-linear empirical erosion laws results in a very strong localization of rapid glacier erosion in regions of fast flow.

• This high spatio-temporal variability hints at a complex relationship between climate and glacier erosion, so a highly variable response of glacier erosion to climate should be expected.

• Increased gravitational drag due to surface profile steepening provides a mechanism for accelerated erosion during deglaciation periods, irrespective of surface meltwater penetration to the glacier bed.

• If a non-linear glacier erosion law is used, the climate-induced slowdown of erosion counterbalances glacier expansion, so Alpine-wide glacier erosion volumes do not correlate with the ice volume.

• Rapid glacier erosion is restricted to low elevation during stages of glacier advance and maximum glaciation but propagates up-valley to higher elevations during periods of glacier retreat.

• The diversity of the Alpine glacier erosion landscape appears to be the time-transgressive signature of a variety of glacial stages ranging from pan-Alpine ice sheets to modern-style mountain glaciers.