Articles | Volume 72, issue 2
Research article
21 Aug 2023
Research article |  | 21 Aug 2023

Subglacial hydrology from high-resolution ice-flow simulations of the Rhine Glacier during the Last Glacial Maximum: a proxy for glacial erosion

Denis Cohen, Guillaume Jouvet, Thomas Zwinger, Angela Landgraf, and Urs H. Fischer
How to cite

Cohen, D., Jouvet, G., Zwinger, T., Landgraf, A., and Fischer, U. H.: Subglacial hydrology from high-resolution ice-flow simulations of the Rhine Glacier during the Last Glacial Maximum: a proxy for glacial erosion, E&G Quaternary Sci. J., 72, 189–201,, 2023.


At the Last Glacial Maximum (LGM), the Rhine Glacier complex (Rhine and Linth glaciers) formed large piedmont lobes extending north into the Swiss and German Alpine forelands. Numerous overdeepened valleys there were formed by repeated glaciations. A characteristic of these overdeepened valleys is their location close to the LGM ice margin, away from the Alps. Numerical models of ice flow of the Rhine Glacier indicate a poor fit between the sliding distance, a proxy for glacial erosion, and the location of these overdeepenings. Calculations of the hydraulic potential based on the computed time-dependent ice surface elevations of the Rhine Glacier lobe obtained from a high-resolution thermo-mechanically coupled Stokes flow model are used to estimate the location of subglacial water drainage routes. Results indicate that the subglacial water discharge is high and focused along glacial valleys and overdeepenings when water pressure is equal to the ice overburden pressure. These conditions are necessary for subglacial water to remove basal sediments, expose fresh bedrock, and favor further erosion by quarrying and abrasion. Knowledge of the location of paleo-subglacial water drainage routes may be useful to understand patterns of subglacial erosion beneath paleo-ice masses that do not otherwise relate to the sliding of ice. Comparison of the erosion pattern from subglacial meltwater with those from quarrying and abrasion shows the importance of subglacial water flow in the formation of distal overdeepenings in the Swiss lowlands.


Während des letzteiszeitlichen Maximums (LGM) kam es zur Bildung von grossen Vorlandloben des Rheingletschersystems (Rhein- und Linthgletscher), die sich nordwärts in das Schweizer und deutsche Alpenvorland erstreckten. Durch wiederholte Vergletscherungen wurden dort zahlreiche übertiefte Täler ausgeschürft. Ein Merkmal dieser übertieften Tälern ist deren Lage in der Nähe des LGM-Eisrands weit weg von den Alpen. Jedoch zeigen numerische Modellierungen des Eisfliessens des Rheingletschers eine schlechte Übereinstimmung der Gleitdistanz, Proxyindikator für glaziale Erosion, mit der Lage dieser Übertiefungen. Deshalb werden Berechnungen des hydraulischen Potenzials basierend auf zeitabhängigen Höhen der Eisoberfläche der Rheingletscherlobe, welche von einem hoch aufgelösten thermo-mechanisch gekoppelten Modell für Stoke’sches Fliessen resultieren, benutzt, um die Lage der Entwässerungsrouten unter dem Eis abzuschätzen. Die Resultate deuten darauf hin, dass der subglaziale Wasserabfluss gross ist und entlang glazialen Tälern und Übertiefungen geführt wird, wenn der Wasserdruck dem Eisüberlagerungsdruck entspricht. Dies sind notwendige Bedingungen, unter denen basale Sedimente wegtransportiert werden und frischer Fels freigelegt wird, um weitere glaziale Erosion zu begünstigen. Somit ist die Kenntnis der Lage von subglazialen Paläo-Entwässerungsrouten nützlich, um die Erosionsmuster unter Paläo-Gletschern zu verstehen, die nicht mit der Gleitbewegung des Eises in Verbindung gebracht werden können. Ein Vergleich der durch subglaziale Schmelzwässer erzeugten Erosionsmuster mit jenen, die durch direkte Gletschererosion entstanden sind, zeigt die Wichtigkeit der subglazialen Schmelzwasserflüsse für die Entstehung von Übertiefungen im alpenfernen Schweizer Vorland auf.

1 Introduction

At the Last Glacial Maximum (LGM), the Rhine Glacier complex (Cohen et al.2018), combining the Rhine and Linth glaciers, descended well into the Swiss and German Alpine forelands, forming two large piedmont lobes. The Rhine Glacier lobe covered present-day Lake Constance and land to the north in southern Germany, and the Linth Glacier lobe advanced beyond the city of Zurich, Switzerland. Repeated glaciations since the Middle–Late Pleistocene (Preusser et al.2011; Ellwanger et al.2011) have carved numerous landforms in the Alps such as deep alpine valleys with prominent horns and ridges. In the forelands, the passage of glaciers left numerous imprints on the landscape such as terminal moraines, outwash deposits, and erratics (e.g., Schlüchter1988, 2004; Keller and Krayss2005a; Preusser et al.2007; Beckenbach et al.2014; Gaar et al.2019); drumlin fields (Kamleitner2022); tunnel valleys (Reber and Schlunegger2016); and overdeepened valleys now filled with sediments such as the Thur, the Glatt, and the Aare valleys (e.g., Preusser et al.2011; Dehnert et al.2012; Dürst Stucki and Schlunegger2013) or with water such as Lake Constance and Lake Zurich. A surprising characteristic of these overdeepenings along these valleys is their proximity to the terminal position of past glacial maxima in a region where one would expect that the erosive power of ice, often associated with the rate of sliding at the base (e.g., Hallet1981; Humphrey and Raymond1994; MacGregor et al.2000; Koppes et al.2015; Herman et al.2015), would be smaller in comparison to that of large alpine valleys where ice fluxes and sliding speeds are higher. The existence of these overdeepenings in a distal-foreland setting thus remains puzzling.

Previous models of ice flow of the Rhine Glacier (Cohen and Jouvet2017; Haeberli et al.2020; Fischer et al.2021; Seguinot and Delaney2021) indicate that the sliding distance, a proxy for glacial erosion based on the time-integrated sliding speed, or a power of it, during the advance and retreat of the Rhine Glacier lobe into the foreland, does not correlate well with the location of existing overdeepenings (Fig. 1) owing to the short time period over which ice covered these distal forelands. More recent work by Egholm (2022) using a landscape evolution model indicates that the formation of distal overdeepenings is possible if the rock material is significantly weaker there than the material in the Alpine massif. The model of Egholm (2022) used empirical parametrizations of abrasion and quarrying (Ugelvig et al.2018) that were fitted against theoretical models (Hallet1981; Iverson2012). Egholm (2022)'s results, however, were mostly obtained using a steady-state climate model in which ice occupation time over the Swiss Alpine foreland was artificially long, in contrast with the geomorphic record for the Rhine Glacier (e.g., Keller and Krayss2005b) that suggests a maximum extent lasting a couple of thousand of years. This dilemma begs the question as to what other subglacial glaciological process could have contributed to the erosion of the northern Alpine forelands and the formation of these distal overdeepenings.

Several authors have noted the importance of subglacial water in glacial erosion (see review by Alley et al.2019). Erosion by quarrying (Hooke1981; Iverson1991; Alley et al.1999), whereby blocks of rock are plucked from the glacier bed by the moving ice, depends on subglacial water pressure fluctuations. Field experiments (Cohen et al.2006) and models of quarrying (Iverson1991, 2012; Hallet1996; Beaud et al.2014; Ugelvig et al.2018) have shown a link between fluctuations in water pressure and the rate of opening of pre-existing bedrock cracks and fractures. Furthermore, large-scale investigations of the erosive power of glaciers (Koppes et al.2015) indicate that temperate basal conditions, which imply the presence of subglacial water, are indicative of higher erosion rates, particularly at middle to high latitudes where seasonal and diurnal temperature variations are high and thus increase the possibility of significant water pressure fluctuations. Subglacial water is also key to the removal of sediments eroded and deposited at the base of glaciers (e.g., Hallet1996) that would otherwise blanket the bedrock, preventing any erosion process from occurring (Alley et al.2019). Indeed, freeze-on conditions along the adverse slope of overdeepenings inhibit subglacial water flow and sediment removal, leading to decrease in erosion there (e.g., Hooke1991; Alley et al.2003; Werder2016). Thus, subglacial water plays an important role in the erosive power of glaciers.

Lacking a clear link between sliding ice (or motion) and observed locations of overdeepenings in the Swiss Alpine forelands underneath the paleo-Rhine Glacier and paleo-Linth Glacier lobes, we seek to find whether subglacial water drainage distribution under the paleo-Rhine Glacier complex could help us understand the location of these overdeepenings. Our approach uses the numerical ice-flow model described in Cohen et al. (2018) but with a new climate signal to drive the ice model and a finer horizontal discretization to better resolve the ice flow. Calculated ice surface elevation, together with basal topography, allows us to compute the hydraulic potential and hydraulic head and the associated flow accumulation area during the course of advance and retreat of the Rhine Glacier around the LGM. Time integration of the flow accumulation area yields a metric for glacial erosion by subglacial water to test against the positions of observed overdeepenings. The underlying assumption in this analysis is that larger values of the time-integrated flow accumulation area indicate zones with significant subglacial water flow, which, over seasonal and annual timescales, have high potential for subglacial water fluxes that would enhance sediment removal and increase glacial erosion. We then compare patterns of relative erosion by quarrying, abrasion, and subglacial water flow to show the importance of subglacial water in the formation of distal overdeepenings in the Swiss lowlands.

The remainder of this paper is divided into three parts. Section 2 presents the methodology to compute the time-integrated flow accumulation area and the erosion by quarrying, abrasion, and subglacial water; Sect. 3 discusses the source of the data used for the computation; finally, Sect. 4 reports results and discusses them.

Figure 1Map of (a) sliding speed and (b) sliding distance calculated by integrating the sliding speed over the advance and the retreat of the Rhine Glacier system during the last glacial cycle (see Cohen2017; Cohen and Jouvet2017; and Cohen et al.2018, for details). The location of overdeepenings is indicated in orange. Reproduced from Fischer et al. (2021).

2 Methods

Subglacial water flow paths are calculated using the method described in Chu et al. (2016) (see also Arnold et al.1998; Willis et al.2012; Livingstone et al.2013) and used in several studies to estimate subglacial water drainage locations in Greenland and Antarctica (e.g., Le Brocq et al.2009; see also references in Chu et al.2016) and in paleo-ice masses such as the Fennoscandian Ice Sheet (e.g., Shackleton et al.2018). More sophisticated models of subglacial hydrology (e.g., GlaDS; Werder et al.2013) coupled to ice-flow models (e.g., Gagliardini and Werder2018, with Elmer/Ice) are too computationally expensive to be applied to the long timescales and the large complex topography of the paleo-Rhine Glacier system.

Water flows down the hydraulic potential gradient, and the water flux, Q, can be calculated as (Shreve1972)

(1) Q = - k ϕ ,

where k is the hydraulic conductivity of the subglacial water system and ϕ is the hydraulic potential, defined as

(2) ϕ = ρ w g Z b + p w .

Here ρw=1000 kg m−3 is the density of water (constant), g=9.81 m s−2 is the acceleration due to gravity, Zb is the elevation of the basal surface (glacier bed topography), and pw is the water pressure. In glaciological applications, the water pressure at the bottom of the ice is often not known. In paleo-glaciological applications, no subglacial water pressure data exist, so it is often not computed (e.g., Cohen et al.2018). Here we assume that the water pressure at the bed of the glacier is a fraction of the ice overburden pressure; i.e.,

(3) p w = f p i ,

where f is known as the flotation factor and

(4) p i = ρ i g H ,

with ρi=917 kg m−3 being the ice density (constant) and H the ice thickness (variable). Substituting Eqs. (3) and (4) into Eq. (2) yields

(5) ϕ = ρ w g Z b + f ρ i g H ,

where the hydraulic potential now only depends on the basal surface topography and the ice thickness. As an alternative we can use the hydraulic head, h, defined as

(6) h = Z b + f ρ i ρ w H .

The value of the flotation factor f depends on the basal water pressure. When water pressure equals ice overburden pressure, f=1. Values greater than 1 are observed during times of intense surface melt that bring water through crevasses and moulins directly to the bed of the glacier (Meierbachtol et al.2013; Wright et al.2016). Since our ice-flow model (Cohen et al.2018) does not simulate subglacial water pressures and given the uncertainties in both the basal surface during the last glacial cycle and estimates of the ice surface, we assume spatially fixed values of f over the entire glacier system and compute the hydraulic potential for three cases – f=0.6, f=1, and f=1.1 – to study the effects of water pressure on subglacial water drainage paths. Low values of f imply that basal topography dominates, while higher values emphasize the ice surface topography (see Eq. 6).

The subglacial water drainage paths can be obtained by computing the upstream accumulation area at each of the cells of the raster of the hydraulic head (Flowers and Clarke1999; Fischer et al.2005; Chu et al.2016; Pitcher et al.2016) using the D algorithm of Tarboton (1997) to compute flow directions. Water drainage flow paths are continuous paths with high values of upstream accumulation areas.

On the long timescales over which erosion occurs (on the order of thousands of years), the ice thickness (H in Eq. 6) evolves with advances and retreats of the glacier. To include these transient effects, we compute the subglacial water flow paths at each available time step for the ice surface (every 10 years) and integrate the computed upstream accumulation area over time to obtain a time-integrated view of subglacial water drainage paths.

Estimates of erosion by subglacial water flow are calculated using the method of Kirkham et al. (2022) based on a model of subglacial channel erosion by Carter et al. (2017). Here we only compute the erosive component (see Kirkham et al.2022, Eq. 3),

(7) E ˙ s = K 1 ν s α cc max τ cc - τ k , 0 g ρ s - ρ w D 15 3 / 2 ,



Here τcc is the channel shear stress, τk is the critical shear stress, νs is the mean sediment velocity, uc is the water velocity in the channel, fcc=0.07 m-2/3 s−2 denotes the channel roughness parameters, ρs=2600 kg m−3 is the solid density of sediment grains, D15=0.1 mm is the 15th-percentile grain size, and μw=1.78×10-3 Pa s is the viscosity of water. The water flow velocity in the channel, uc, is computed from the flow accumulation area weighted by the available meltwater at each cell, which is estimated by summing the surface melt rate and the basal melt rate. The meltwater is then scaled by multiplying it by the grid cell area divided by the cross-sectional area of the channel, assumed to be 30 m wide and 2 m deep, to obtain the channel water velocity. Although changing the channel size would affect the relationship between τcc and τk and thus where erosion occurs, it does not affect the overall results of erosion patterns as our calculations are meant to only look at erosion patterns and not calculate erosion magnitudes. The surface melt rate is calculated in the ablation area of the glacier and is proportional to the number of positive degree days (PDDs) (Hock2003). It is assumed to be zero in the accumulation area. The basal melt rate is computed from the geothermal heat flux map of Medici and Rybach (1995). Frictional heat generated by ice sliding over the substrate is an order of magnitude smaller than either surface or basal melt and is thus neglected.

To compare patterns of subglacial erosion due to ice motion, we also estimate erosion by quarrying, E˙q, and abrasion, E˙a, given, respectively, by (see Ugelvig et al.2016, 2018)


where pe=pw-pi is the effective pressure, vs is the glacier sliding speed, and Kq and Ka are constants.

All erosion rates are time-integrated to yield a total erosion over the course of the advance and retreat of the glacier during the LGM. Because of uncertainties in model and parameter values, a direct comparison between quarrying, abrasion, and subglacial meltwater erosion is not possible. Instead, these quantities are used to look at relative patterns of erosion by normalizing them. This makes the exact values of the parameters K1 (Eq. 7) and Kq and Ka (Eqs. 11 and 12, respectively) irrelevant.

3 Data

Only two surfaces are required to compute the hydraulic potential and subglacial water flow paths: the basal surface (Zb in Eq. 6) and the ice surface from which the ice thickness, H, can be obtained (see Eq. 6). The ice surface is computed from ice-flow simulations of the paleo-Rhine Glacier complex around the LGM using Elmer/Ice, an open-source, finite-element multi-physics Fortran code (Gagliardini et al.2013; Ruokolainen et al.2020). Elmer/Ice computes the transient three-dimensional ice velocities, temperatures, and pressures as well as the elevation of the ice surface based on fully coupled thermo-mechanical equations by solving the non-linear Stokes flow equation together with the heat equation and the mass balance equation at the ice surface (see Cohen et al.2018, for details). The mass balance is based on the positive-degree-day model (Hock2003) that computes ice melt as a function of air temperature. PDDs are computed using the method of Calov and Greve (2005) over month-long intervals using factors of 8 and 3 mm C−1 d−1 for ice and snow, respectively. Air temperature and precipitation are obtained by linear interpolation using a glacial-index (GI) approach (Sutter et al.2019) between LGM and pre-industrial (PI) climate states. The GI follows the δ18O signal of the Antarctica EPICA ice core (Jouzel et al.2007), rescaled to between 0 (no ice, PI) and 1 (maximum ice volume at the LGM). The LGM climate state provides 2 km resolution maps of the temperature and precipitation over the entire Alps computed from a regional climate model by dynamical downscaling of a global Earth system model (see details in Velasquez et al.2020, 2022; Russo et al.2022; Buzan et al.2023).

3.1 Basal surface

The basal surface is the present-day topography minus the computed ice thickness of major present-day glaciers using the model of Farinotti et al. (2019). It remains unchanged during the ice-flow simulation. Figure 2 shows the basal topography used in the model simulation.

Figure 2Basal topography with outline of overdeepened basins in black and extent of ice at the LGM with orange circles (Ehlers and Gibbard2008).​​​​​​​

3.2 Ice surface

The ice surface evolves with time during the advance of the Rhine Glacier towards its maximum LGM position and during its retreat. Figure 3 shows three ice surface elevations at 27 ka, at 22.9 ka (LGM), and at 14.2 ka. Ice surface elevations are available every 10 years from 27 to 14.2 ka.

Figure 3Ice surface elevation in blue at (a) 27 ka (start of simulation), (b) 22.9 ka (LGM), and (c) 14.2 ka (end of simulation). LGM ice extent is shown with orange circles (Ehlers and Gibbard2008).

3.3 Calculation of hydraulic potential and flow accumulation area

The hydraulic potential is first calculated on the unstructured triangular mesh of the Elmer/Ice finite-element model at every time step (every 10 years). Resolution in the mesh ranges from 800 to 1100 m, with increased resolution in areas of steeper slope gradients. The nodal values of the hydraulic potential are then linearly interpolated to a uniform raster at a resolution of 500 m. The calculation of the flow accumulation area is performed on the hydraulic potential using the open toolbox SAGA (System for Automated Geoscientific Analyses; Conrad et al.2015) and summed over time to compute the time-integrated values.

4 Results and discussion

Figure 4 shows the subglacial water drainage routes based on the flow accumulation area of the hydraulic potential (or hydraulic head) at three different times for flotation values of f=0.6 (Fig. 4, left panels), f=1.0 (Fig. 4, center panels), and f=1.1 (Fig. 4, right panels).

Figure 4The log 10 of computed flow accumulation area (in m2) in blue at three different times – (a, b, c) 27 ka, (d, e, f) 22.9 ka (LGM), and (g, h, i) 14.2 ka – based on the hydraulic potential for a flotation factor of (a, d, g) f=0.6, (b, e, h) f=1.0, and (c, f, i) f=1.1. Hillshade topography and ice thickness in red shades at the times indicated are shown in the background. Actual LGM ice extent is shown with orange circles. Modeled ice extent is shown with a brown line. Black outlines indicate locations of overdeepenings.

The flow accumulation area can be understood as a proxy for subglacial water flux that passes through a specific drainage channel because water flux to the glacier is an increasing function of the upstream flow accumulation area.

When f=0.6 (Fig. 4, left panels), subglacial water drainages are mostly concentrated in the Alpine Rhein Valley in the Alps, in Lake Constance, and along major valleys in the forelands such as the drainage of Lake Zurich. Drainages are constrained by basal topography, following valley morphology. Lake Constance (a zone of constant elevation), where the Alpine Rhein Glacier flows, concentrates a significant flux of subglacial water according to the model. Water exits the Rhine Glacier lobe to the north via existing valleys in the German Alpine foreland, following the Überlingen branch of present-day Lake Constance. In the Swiss Alpine foreland, water drains out of Lake Constance via the Untersee and also by connecting to the Thur Valley, which coincides with an overdeepening there. In the Linth Glacier lobe, subglacial water also flows along Lake Zurich and the Limmat and Glatt valleys, where two overdeepenings are located (Preusser et al.2010; Buechi et al.2018). The Reuss and Aare overdeepenings are also pathways for subglacial water drainages further to the west. The drainage pattern does not change significantly with time (see Fig. 4a, c, e) except at t=14.2 ka when the ice has retreated halfway out of Lake Constance and completely out of the distal overdeepenings in the Swiss forelands.

The drainage pattern is significantly different when f=1.0 (Fig. 4, center panels), equivalent to zero effective pressure at the base. Water now drains out of the Lake Constance basin to the north in many concentrated flow channels, such as the Schussen Valley, against an adverse bedrock slope (see basal topography in Fig. 2) and where a few overdeepenings have been noted (Ellwanger et al.2011). Subglacial flow is also more important along many valleys in the foreland, both in the Rhine Glacier (e.g., Thur) and in the Linth Glacier (e.g., Reuss) lobes. Because of higher water pressure, the hydraulic potential gradient conforms less to the basal (or bed) topography, and subglacial water moves more along the gradient determined by the slope of the glacier surface; subglacial water paths no longer follow topographic lows.

At a yet higher flotation value (f=1.1, Fig. 4, right panels), when water pressure exceeds ice overburden pressure, subglacial drainages out of the Alpine Rhine are also diverted to the north across Lake Constance and then ascending the bedrock slope north of it. Water drainage, however, occurs over many more subglacial channels that fan out radially from the Lake Constance basin. Because of the high water pressures, drainages are no longer constrained by basal topography and can move from one valley to the next over topographic high points such as from the Thur to the Lower Rhine overdeepenings out of the Untersee. High flotation values cause water drainage paths to cut across valleys and overdeepenings, such as, for example, in the Thur Valley, in the Rhine Glacier lobe, and in the Reuss and Linth lobes. In the Alps, water drainage paths are no longer constrained to the valley centers but meander on the sides of the valley and along the valley walls (contrast Fig. 4d and f), even cutting across mountain ranges during retreat (see Fig. 4i).

Basal conditions where the value of flotation is equal to 1 or greater are unlikely to occur over the entire glacier basal surface and for significant lengths of time. High basal water pressure arising from increased surface melting would promote the development of an efficient drainage system that would tend to reduce basal water pressure (e.g., Bartholomaus et al.2008; Sundal et al.2011; Andrews et al.2014). These basal conditions, despite their limited duration, could have the potential to produce significant changes in the basal geomorphology, for example, because of substantial increases in the capacity of subglacial water streams to transport sediments and in their ability to further erode bedrock by processes of abrasion and quarrying associated with both larger fluctuations in basal water pressure and increased basal sliding speed (e.g., van de Wal et al.2008; Schoof2010; Bartholomew et al.2010; Hewitt2013).

The time integration of subglacial water flow paths over the 12.8 kyr of the simulation covering advance and retreat of the Rhine Glacier around the LGM is shown in Fig. 5 for the three values of flotation: f=0.6, 1, and 1.1.

Figure 5The log 10 of the time-integrated flow accumulation area (in m2) over 12.8 kyr (in blue) for (a) f=0.6, (b) f=1, and f=1.1 with hillshade topography and ice thickness at the LGM (red colors) in the background. LGM ice extent is shown with orange circles. Model ice extent is shown with a brown line. Black outlines indicate location of overdeepenings.

Clear differences exist between a moderate flotation value (f=0.6) and flotation values with zero effective pressure (f=1) or where water pressure exceeds ice overburden pressure (f=1.1). Although high water pressures in excess of hydrostatic pressure have been noted in both Alpine (e.g., Flowers and Clarke1999) and Greenland (e.g., Banwell et al.2013; Lindbäck et al.2015) settings and may be occurring in overdeepenings (Lawson et al.1998), they may not be representative of long-timescale average values in large areas of the Rhine Glacier lobe. Subglacial water flow drainages with high water pressure bypass some topographic relief, producing near-straight channels that cut across relief, such as in the northern part of the Rhine Glacier lobe. These high-water-pressure subglacial drainage paths also cut across overdeepenings in the Thur, Glatt, and Limmat overdeepenings in the Rhine Glacier and Linth Glacier lobes. A value of f equal to unity yields the closest match with the actual mapping of overdeepenings, although many overdeepenings are also well matched for f=0.6, particularly in the Linth Glacier lobe and the lower Thur Valley. Overdeepened regions north of Lake Constance require higher flotation values to be occupied by a subglacial water drainage path. When water pressure exceeds ice overburden pressure such as when f=1.1, subglacial drainages are too distributed and no longer only conform to valleys, a situation that may occur during periods of significant melt (such as during retreat) but is likely not representative of the situation over longer timescales.

High values of the time-integrated flow accumulation area represent drainage pathways with large fluxes of water. Many of these high-discharge zones correspond to several locations of overdeepenings such as the Schussen, Thur, Limmat, Reuss, and Aare overdeepened valleys. In the Aare Valley, a tunnel valley system has been identified beneath the city of Bern (Dürst Stucki and Schlunegger2013; Reber and Schlunegger2016) with a branching pattern that resembles the pattern observed in our simulation. The presence of this tunnel valley system has been associated with an area of high water pressure. Our simulation also indicates that a flotation value of 1 is necessary for significant subglacial water drainage to occur in this area.

Figure 6 shows patterns of erosion for quarrying, abrasion, and subglacial water flow. All values of erosion have been scaled to between 0 and 1 by dividing the computed values by their respective maximum value. All erosion types have been time-integrated over the course of the advance and retreat of the ice during the LGM.

Figure 6Normalized erosion by (a) quarrying, (b) abrasion, and (c) subglacial meltwater flow. Background shows hillshade topography. LGM margin is shown with orange circles; modeled maximum LGM extent is shown with a brown curve. Black curves indicate areas of overdeepenings. The red curve in (c) is the equilibrium line at the LGM separating accumulation (south) from ablation (north) in the Alps.

A clear difference in the erosion pattern emerges between erosion due to ice motion (quarrying and abrasion) and erosion due to subglacial meltwater flow. Erosion by ice motion is greatest at the base of alpine valleys where sliding is fastest and where the ice occupation time is longest. Limited durations of ice cover and lower sliding speeds near the margin limit erosion by ice motion there. Erosion by quarrying (Fig. 6a) shows more variability in the Rhine, Linth, and Reuss lobes owing to its cubic dependence on the ice thickness via the effective pressure term (see Eq. 11). Any variation in ice thickness along or across overdeepenings amplifies erosion by quarrying, resulting in a significant change in the total amount of quarrying in these areas. In contrast, erosion by abrasion (Fig. 6b), which only depends on the square of the sliding speed (see Eq. 11), is more uniform in these lobes. Erosion by subglacial meltwater flow (Fig. 6c) is significantly different: it is almost entirely focused near the ice margin, particularly to the northwest in the Reuss and Linth lobes, the Thur Valley, and the Argen lobe (see Fig. 2). This is due to the large available surface melt in these areas that have larger ablation areas. This is visible by comparing, for example, the position of the equilibrium line (Fig. 6c, red line) to the LGM margin (Fig. 6c, brown line). In the northwestern lobes and in the northeastern part of the Rhine lobe, the area for available surface meltwater is significantly larger then in the main Rhine lobe, explaining why most erosion by subglacial meltwater occurs preferentially there. Also, channel shear stress is more likely to exceed critical shear stress in areas where subglacial channel velocity is higher (see uc in Eq. 8), which occurs where more surface melt is available and thus increases the potential for erosion by subglacial meltwater.

The different erosion patterns due to ice motion and subglacial meltwater suggest that the existence of distal-foreland overdeepenings in the Swiss lowlands is driven, at least in part, by erosion by subglacial meltwater. Furthermore, high water flux during high-water-pressure events along these valleys and the associated fluctuations in water pressure in the channels could have increased the rate of glacial erosion by quarrying (e.g., Cohen et al.2006), a dependence that is not included in the simple quarrying erosion law in Eq. (11). This may have been further enhanced by pre-existing rock weaknesses due, for example, to a higher density of joints or faulting (e.g., Dühnforth et al.2010; Hooyer et al.2012) such as is found in the Lake Constance area (Fabbri et al.2021). The patterns of subglacial water drainage paths computed from numerical simulations of ice flow and ice surface elevation indicate a possible link between the subglacial water flux and the existence of overdeepenings. The location of these overdeepenings suggests that these overdeepenings are controlled, in part, by higher rock erodibility in these areas, which, during glacial times, were preferentially excavated by the combined effects of water and ice.

5 Conclusions

We present a model of subglacial water drainages underneath the Rhine Glacier during the advance and retreat around the LGM. The model is based on high-resolution ice surfaces obtained from a numerical ice-flow model. The model computes ice velocities and ice surface elevations using a thermo-mechanical Stokes flow model coupled to a high-resolution regional climate model to give the most plausible representation of the Rhine Glacier during its advance and retreat around the LGM. The hydraulic potential (or head) calculated from transient elevations of the ice surface together with the basal topography yields a map of subglacial water flux during the LGM. Results of the calculation show that subglacial water drainage is strongly constrained by topography for low values of water pressure (small flotation factor), focusing water flow in valleys in the foreland. When water pressure equals ice overburden pressure (flotation value equal to 1), the subglacial drainages coincide with regions of intense glacial erosion as observed by numerous overdeepenings there. When water pressure exceeds ice overburden pressure and flotation is greater than 1, the pattern of subglacial drainages changes to straighter channels that bypass valleys, mountain ranges, and overdeepened basins. This suggests that the maximum potential for glacial erosion by subglacial water channels is obtained when water pressure underneath the Rhine Glacier and Linth Glacier lobes were, on average, close to the ice overburden pressure. These conditions favor high water flow volumes and maximize the potential for subglacial water fluctuations needed for enhanced erosion by ice and water. Finally, comparison of the erosion pattern from subglacial meltwater with those from quarrying and abrasion shows the importance of subglacial water channels in the formation of distal overdeepenings in the Swiss lowlands.

Code availability

The software package Elmer FEM version 9 which contains Elmer/Ice can be downloaded via Zenodo, (Ruokolainen et al.2023).

Data availability

Data are not publicly available as part of the climate data were obtained from another group (references cited in text) and they have not been made available publicly as of now. Other data sources are not publicly available as they were obtained from Nagra.

Author contributions

The research was conceived by DC, UHF, and AL. The ice-flow modeling was performed by DC with help from TZ and GJ. DC carried out hydrological modeling and interpretation. DC wrote the paper with inputs from all co-authors.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Subglacial erosional landforms and their relevance for the long-term safety of a radioactive waste repository”. It is the result of a virtual workshop held in December 2021.

Financial support

This research has been supported by the National Cooperative for the Disposal of Radioactive Waste (Nagra, grant no. 20632).

Review statement

This paper was edited by Sonja Breuer and reviewed by Meike Bagge and Anders Damsgaard.


Alley, R. B., Strasser, J. C., Lawson, D. E., Evenson, E. B., and Larson, G. J.: Glaciological and geological implications of basal-ice accretion in overdeepenings, Geol. S. Am. S., 337, 1–9,, 1999. a

Alley, R. B., Lawson, D. E., Evenson, E. B., and Larson, G. J.: Sediment, glaciohydraulic supercooling, and fast glacier flow, Ann. Glaciol., 36, 135–141,, 2003. a

Alley, R. B., Cuffey, K. M., and Zoet, L. K.: Glacial erosion: status and outlook, Ann. Glaciol., 60, 1–13,, 2019. a, b

Andrews, L. C., Catania, G. A., Hoffman, M. J., Gulley, J. D., Lüthi, M. P., Ryser, C., Hawley, R. L., and Neumann, T. A.: Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet, Nature, 514, 80–83,, 2014. a

Arnold, N., Richards, K., Willis, I., and Sharp, M.: Initial results from a distributed, physically based model of glacier hydrology, Hydrol. Process., 12, 191–219,<191::AID-HYP571>3.0.CO;2-C, 1998. a

Banwell, A. F., Willis, I. C., and Arnold, N. S.: Modeling subglacial water routing at Paakitsoq, W Greenland, J. Geophys. Res.-Earth Surf., 118, 1282–1295,, 2013. a

Bartholomaus, T. C., Anderson, R. S., and Anderson, S. P.: Response of glacier basal motion to transient water storage, Nat. Geosci., 1, 33–37,, 2008. a

Bartholomew, I., Nienow, P., Mair, D., Hubbard, A., King, M. A., and Sole, A.: Seasonal evolution of subglacial drainage and acceleration in a Greenland outlet glacier, Nat. Geosci., 3, 408–411,, 2010. a

Beaud, F., Flowers, G. E., and Pimentel, S.: Seasonal-scale abrasion and quarrying patterns from a two-dimensional ice-flow model coupled to distributed and channelized subglacial drainage, Geomorphology, 219, 176–191,, 2014. a

Beckenbach, E., Müller, T., Seyfried, H., and Simon, T.: Potential of a high-resolution DTM with large spatial coverage for visualization, identification and interpretation of young (Würmiam) glacial geomorphology – a case study from Oberschwaben (southern Germany), Quaternary Science Journal, 63, 107–129,, 2014. a

Buechi, M. W., Graf, H. R., Haldimann, P., Lowick, S. E., and Anselmetti, F. S.: Multiple Quaternary erosion and infill cycles in overdeepened basins of the northern Alpine foreland, Swiss J. Geosci., 111, 133–167,, 2018. a

Buzan, J. R., Russo, E., Kim, W. M., and Raible, C. C.: Winter sensitivity of glacial states to orbits and ice sheet heights in CESM1.2, EGUsphere [preprint],, 2023. a

Calov, R. and Greve, R.: Correspondence: A semi-analytical solution for the positive degree-day model with stochastic temperature variations, J. Glaciol, 51, 173–175, 2005. a

Carter, S. P., Fricker, H. A., and Siegfried, M. R.: Antarctic subglacial lakes drain through sediment-floored canals: theory and model testing on real and idealized domains, The Cryosphere, 11, 381–405,, 2017. a

Chu, W., Creyts, T. T., and Bell, R. E.: Rerouting of subglacial water flow between neighboring glaciers in West Greenland, J. Geophys. Res.-Earth Surf., 121, 925–938,, 2016. a, b, c

Cohen, D.: Numerical Reconstruction of the Rhine Glacier at the Last Glacial Maximum, Tech. rep., Nagra Arbeitsbericht NAB 17–25,, 2017. a

Cohen, D. and Jouvet, G.: Transient Simulations of the Rhine Glacier over the Last Glacial Cycle, Tech. rep., Nagra Arbeitsbericht NAB 17–47, 2017. a, b

Cohen, D., Hooyer, T. S., Iverson, N. R., Thomason, J. F., and Jackson, M.: Role of transient water pressure in quarrying: A subglacial experiment using acoustic emissions, J. Geophys. Res.-Earth Surf., 111, F03006,, 2006. a, b

Cohen, D., Gillet-Chaulet, F., Haeberli, W., Machguth, H., and Fischer, U. H.: Numerical reconstructions of the flow and basal conditions of the Rhine glacier, European Central Alps, at the Last Glacial Maximum, The Cryosphere, 12, 2515–2544,, 2018. a, b, c, d, e, f

Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L., Wehberg, J., Wichmann, V., and Böhner, J.: System for Automated Geoscientific Analyses (SAGA) v. 2.1.4, Geosci. Model Dev., 8, 1991–2007,, 2015. a

Dehnert, A., Lowick, S. E., Preusser, F., Anselmetti, F. S., Drescher-Schneider, R., Graf, H. R., Heller, F., Horstmeyer, H., Kemna, H. A., Nowaczyk, N. R., Züger, A., and Furrer, H.: Evolution of an overdeepened trough in the northern Alpine Foreland at Niederweningen, Switzerland, Quaternary Sci. Rev., 34, 127–145,, 2012. a

Dühnforth, M., Anderson, R. S., Ward, D., and Stock, G. M.: Bedrock fracture control of glacial erosion processes and rates, Geology, 38, 423–426,, 2010. a

Dürst Stucki, M. and Schlunegger, F.: Identification of erosional mechanisms during past glaciations based on a bedrock surface model of the central European Alps, Earth Planet. Sci. Lett., 384, 57–70, 2013. a, b

Egholm, D.: Controls on Overdeepening Formation in a Distal Foreland Setting, Tech. rep., Nagra Arbeitsbericht NAB 22-17, 2022. a, b, c

Ehlers, J. and Gibbard, P.: Extent and chronology of Quaternary glaciation, Episodes, 31, 211–218, 2008. a, b

Ellwanger, D., Wielandt-Schuster, U., Franz, M., and Simon, T.: The Quaternary of the southwest German Alpine Foreland (Bodensee–Oberschwaben, Baden-Württemberg, Southwest Germany), Quaternary Sci. J., 60, 306–328, 2011. a, b

Fabbri, S., Affentranger, C., Krastel, S., Lindhorst, K., Wessels, M., Madritsch, H., Allenbach, R., Herwegh, M., Heuberger, S., Wielandt-Schuster, U., Pomella, H., Schwestermann, T., and Anselmetti, F.: Active Faulting in Lake Constance (Austria, Germany, Switzerland) Unraveled by Multi-Vintage Reflection Seismic Data, Front. Earth Sci., 9,, 2021. a

Farinotti, D., Huss, M., Fürst, J. J., Landmann, J., Machguth, H., Maussion, F., and Pandit, A.: A consensus estimate for the ice thickness distribution of all glaciers on Earth, Nat. Geosci., 12, 168–173,, 2019. a

Fischer, U. H., Braun, A., Bauder, A., and Flowers, G. E.: Changes in geometry and subglacial drainage derived from digital elevation models: Unteraargletscher, Switzerland, 1927–97, Ann. Glaciol., 40, 20–24,, 2005. a

Fischer, U. H., Bebiolka, A., Brandefelt, J., Cohen, D., Harper, J., Hirschorn, S., Jensen, M., Kennell, L., Liakka, J., Näslund, J.-O., Normani, S., Stück, H., and Weitkamp, A.: Chapter 11 – Radioactive waste under conditions of future ice ages, in: Snow and Ice-Related Hazards, Risks, and Disasters (Second Edition), edited by: Haeberli, W. and Whiteman, C., Hazards and Disasters Series, 323–375, Elsevier, 2nd Edn.,, 2021. a, b

Flowers, G. E. and Clarke, G. K. C.: Surface and bed topography of Trapridge Glacier, Yukon Territory, Canada: digital elevation models and derived hydraulic geometry, J. Glaciol., 45, 165–174,, 1999. a, b

Gaar, D., Graf, H. R., and Preusser, F.: New chronological constraints on the timing of Late Pleistocene glacier advances in northern Switzerland, E&G Quaternary Sci. J., 68, 53–73,, 2019. a

Gagliardini, O. and Werder, M. A.: Influence of increasing surface melt over decadal timescales on land-terminating Greenland-type outlet glaciers, J. Glaciol., 64, 700–710,, 2018. a

Gagliardini, O., Zwinger, T., Gillet-Chaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a new-generation ice sheet model, Geosci. Model Dev., 6, 1299–1318,, 2013. a

Haeberli, W., Fischer, U. H., Cohen, D., and Schnellmann, M.: Radioaktive Abfälle und Eiszeiten in der Schweiz: Können Gletscher und Permafrost zukünftiger Eiszeiten die langfristige Sicherheit der geplanten Lager beeinflussen?, Wasser Energie Luft, 112, 261–269, 2020. a

Hallet, B.: Glacial Abrasion and Sliding: their Dependence on the Debris Concentration in Basal Ice, Ann. Glaciol., 2, 23–28,, 1981. a, b

Hallet, B.: Glacial quarrying: a simple theoretical model, Ann. Glaciol., 22, 1–8,, 1996. a, b

Herman, F., Beyssac, O., Brughelli, M., Lane, S. N., Leprince, S., Adatte, T., Lin, J. Y. Y., Avouac, J.-P., and Cox, S. C.: Erosion by an Alpine glacier, Science, 350, 193–195,, 2015. a

Hewitt, I. J.: Seasonal changes in ice sheet motion due to melt water lubrication, Earth Planet. Sci. Lett., 371–372, 16–25,, 2013. a

Hock, R.: Temperature index melt modelling in mountain areas, J. Hydrol., 282, 104–115, 2003. a, b

Hooke, R.: Flow law for polycrystalline ice in glaciers: comparison of theoretical predictions, laboratory data, and field measurements, Rev. Geophys. Space. Phys., 19, 664–672, 1981. a

Hooke, R. L.: Positive feedbacks associated with erosion of glacial cirques and overdeepenings, GSA Bulletin, 103, 1104,<1104:PFAWEO>2.3.CO;2, 1991. a

Hooyer, T. S., Cohen, D., and Iverson, N. R.: Control of glacial quarrying by bedrock joints, Geomorphology, 153, 91–101,, 2012. a

Humphrey, N. F. and Raymond, C. F.: Hydrology, erosion and sediment production in a surging glacier: Variegated Glacier, Alaska, 1982–83, J. Glaciol., 40, 539–552,, 1994. a

Iverson, N. R.: Morphology of glacial striae: Implications for abrasion of glacier beds and fault surfaces, GSA Bulletin, 103, 1308–1316,<1308:MOGSIF>2.3.CO;2, 1991. a, b

Iverson, N. R.: A theory of glacial quarrying for landscape evolution models, Geology, 40, 679–682,, 2012. a, b

Jouzel, J., Masson-Delmotte, V., Cattani, O., Dreyfus, G., Falourd, S., Hoffmann, G., Minster, B., Nouet, J., Barnola, J.-M., Chappellaz, J., Fischer, H., Gallet, J. C., Johnsen, S., Leuenberger, M., Loulergue, L., Luethi, D., Oerter, H., Parrenin, F., Raisbeck, G., Raynaud, D., Schilt, A., Schwander, J., Selmo, E., Souchez, R., Spahni, R., Stauffer, B., Steffensen, J. P., Stenni, B., Stocker, T. F., Tison, J. L., Werner, M., and Wolff, E. W.: Orbital and millennial Antarctic climate variability over the past 800 000 years, Science, 317, 793–796, 2007. a

Kamleitner, S.: Reconstructing the evolution and dynamics of Central Alpine glaciers during the Last Glacial Maximum on the basis of their geomorphological footprints and cosmogenic nuclide surface exposure dating, PhD thesis, University of Zürich, 186 pp.,, 2022. a

Keller, O. and Krayss, E.: Der Rhein-Linth-Gletscher im letzten Hochglazial. 1. Teil: Einleitung; Aufbau und Abschmelzen des Rhein-Linth-Gletschers im Oberen Würm, Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich, 150, 19–32, 2005a. a

Keller, O. and Krayss, E.: Der Rhein-Linth-Gletscher im letzten Hochglazial. 2. Teil: Datierung und Modelle der Rhein-Linth-Vergletscherung. Klima-Rekonstruktionen, Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich, 150, 69–85, 2005b. a

Kirkham, J. D., Hogan, K. A., Larter, R. D., Arnold, N. S., Ely, J. C., Clark, C. D., Self, E., Games, K., Huuse, M., Stewart, M. A., Ottesen, D., and Dowdeswell, J. A.: Tunnel valley formation beneath deglaciating mid-latitude ice sheets: Observations and modelling, Quaternary Sci. Rev., 107680,, in press, 2022. a, b

Koppes, M., Hallet, B., Rignot, E., Mouginot, J., Wellner, J. S., and Boldt, K.: Observed latitudinal variations in erosion as a function of glacier dynamics, Nature, 526, 100–103,, 2015. a, b

Lawson, D. E., Strasse r, J. C., Evenson, E. B., Alley, R. B., Larson, G. J., and Arcone, S. A.: Glaciohydraulic supercooling: a freeze-on mechanism to create stratified, debris-rich basal ice: I. Field evidence, J. Glaciol., 44, 547–562,, 1998. a

Le Brocq, A., Payne, A., Siegert, M., and Alley, R.: A subglacial water-flow model for West Antarctica, J. Glaciol., 55, 879–888,, 2009. a

Lindbäck, K., Pettersson, R., Hubbard, A. L., Doyle, S. H., van As, D., Mikkelsen, A. B., and Fitzpatrick, A. A.: Subglacial water drainage, storage, and piracy beneath the Greenland ice sheet, Geophys. Res. Lett., 42, 7606–7614,, 2015. a

Livingstone, S. J., Clark, C. D., Woodward, J., and Kingslake, J.: Potential subglacial lake locations and meltwater drainage pathways beneath the Antarctic and Greenland ice sheets, The Cryosphere, 7, 1721–1740,, 2013. a

MacGregor, K. R., Anderson, R., Anderson, S., and Waddington, E.: Numerical simulations of glacial-valley longitudinal profile evolution, Geology, 28, 1031–1034, 2000. a

Medici, F. and Rybach, L.: Geothermal Map of Switzerland 1995 (Heat Flow Density), Tech. Rep. 30, Swiss Geophysical Commission, 1995. a

Meierbachtol, T., Harper, J., and Humphrey, N.: Basal Drainage System Response to Increasing Surface Melt on the Greenland Ice Sheet, Science, 341, 777–779,, 2013. a

Pitcher, L. H., Smith, L. C., Gleason, C. J., and Yang, K.: CryoSheds: a GIS modeling framework for delineating land-ice watersheds for the Greenland Ice Sheet, GISci. Remote Sens., 53, 707–722,, 2016. a

Preusser, F., Blei, A., Graf, H., and Schlüchter, C.: Luminescence dating of Würmian (Weichselian) proglacial sediments from Switzerland: methodological aspects and stratigraphical conclusions, Boreas, 36, 130–142,, 2007. a

Preusser, F., Reitner, J., and Schlüchter, C.: Distribution, geometry, age and origin of overdeepened valleys and basins in the Alps and their foreland, Swiss J. Geosci., 103, 407–426,, 2010. a

Preusser, F., Graf, H. R., Keller, O., Krayss, E., and Schlüchter, C.: Quaternary glaciation history of northern Switzerland, E&G Quaternary Sci. J., 60, 21,, 2011. a, b

Reber, R. and Schlunegger, F.: Unravelling the moisture sources of the Alpine glaciers using tunnel valleys as constraints, Terra Nova, 28, 202–211,, 2016. a, b

Ruokolainen, J., Malinen, M., Råback, P., Zwinger, T., Pursula, A., and Byckling, M.​​​​​​​: ElmerSolver Manual, Tech. Rep. Online, CSC – IT Center for Science, (last access: 9 August 2023), 2020. a

Ruokolainen, J., Malinen, M., Råback, P., Zwinger, T., Takala, E., Kataja, J., Gillet-Chaulet, F., Ilvonen, S., Gladstone, R., Byckling, M., Chekki, M., Gong, C., Ponomarev, P., van Dongen, E., Robertsen, F., Wheel, I., Cook, S., t7saeki, luzpaz, and Rich_B.​​​​​​​: ElmerCSC/elmerfem: Elmer 9.0 (release-9.0), Zenodo [code],, 2023. a

Russo, E., Fallah, B., Ludwig, P., Karremann, M., and Raible, C. C.: The long-standing dilemma of European summer temperatures at the mid-Holocene and other considerations on learning from the past for the future using a regional climate model, Clim. Past, 18, 895–909,, 2022. a

Schlüchter, C.: A non-classical summary of the Quaternary stratigraphy in the Northern Alpine Foreland of Switzerland, Bulletin de la Société neuchâteloise de Géographie, 32, 143–157, 1988. a

Schlüchter, C.: The Swiss glacial record – A schematic summary, in: Quaternary Glaciations Extent and ChronologyPart I: Europe, edited by: Ehlers, J. and Gibbard, P., Vol. 2, Part 1 of: Developments in Quaternary Sciences, 413–418, Elsevier,, 2004. a

Schoof, C.: Ice sheet acceleration driven by melt supply variability, Nature, 468, 803–806, 2010. a

Seguinot, J. and Delaney, I.: Last-glacial-cycle glacier erosion potential in the Alps, Earth Surf. Dynam., 9, 923–935, 2021. a

Shackleton, C., Patton, H., Hubbard, A., Winsborrow, M., Kingslake, J., Esteves, M., Andreassen, K., and Greenwood, S. L.: Subglacial water storage and drainage beneath the Fennoscandian and Barents Sea ice sheets, Quaternary Sci. Rev., 201, 13–28,, 2018. a

Shreve, R.: Movement of water in glaciers, J. Glaciol, 11, 205–214, 1972. a

Sundal, A. V., Shepherd, A., Nienow, P., Hanna, E., Palmer, S., and Huybrechts, P.: Melt-induced speed-up of Greenland ice sheet offset by efficient subglacial drainage, Nature, 469, 521–524,, 2011. a

Sutter, J., Fischer, H., Grosfeld, K., Karlsson, N. B., Kleiner, T., Van Liefferinge, B., and Eisen, O.: Modelling the Antarctic Ice Sheet across the mid-Pleistocene transition – implications for Oldest Ice, The Cryosphere, 13, 2023–2041,, 2019. a

Tarboton, D. G.: A new method for the determination of flow directions and upslope areas in grid digital elevation models, Water Resour. Res., 33, 309–319,, 1997.  a

Ugelvig, S. V., Egholm, D. L., and Iverson, N. R.: Glacial landscape evolution by subglacial quarrying: A multiscale computational approach, J. Geophys. Res.-Earth Surf., 121, 2042–2068,, 2016. a

Ugelvig, S. V., Egholm, D. L., Anderson, R. S., and Iverson, N. R.: Glacial Erosion Driven by Variations in Meltwater Drainage, J. Geophys. Res.-Earth Surf., 123, 2863–2877,, 2018. a, b, c

van de Wal, R. S. W., Boot, W., van den Broeke, M. R., Smeets, C. J. P. P., Reijmer, C. H., Donker, J. J. A., and Oerlemans, J.: Large and Rapid Melt-Induced Velocity Changes in the Ablation Zone of the Greenland Ice Sheet, Science, 321, 111–113,, 2008. a

Velasquez, P., Messmer, M., and Raible, C. C.: A new bias-correction method for precipitation over complex terrain suitable for different climate states: a case study using WRF (version 3.8.1), Geosci. Model Dev., 13, 5007–5027,, 2020. a

Velasquez, P., Messmer, M., and Raible, C. C.: The role of ice-sheet topography in the Alpine hydro-climate at glacial times, Clim. Past, 18, 1579–1600,, 2022. a

Werder, M., Hewitt, I., Schoof, C., and Flowers, G.: Modeling channelized and distributed subglacial drainage in two dimensions, J Geophys. Res.-Earth Surf., 118, 2140–2158, 2013. a

Werder, M. A.: The hydrology of subglacial overdeepenings: A new supercooling threshold formula, Geophys. Res. Lett., 43, 2045–2052,, 2016. a

Willis, I. C., Fitzsimmons, C. D., Melvold, K., Andreassen, L. M., and Giesen, R. H.: Structure, morphology and water flux of a subglacial drainage system, Midtdalsbreen, Norway, Hydrol. Process., 26, 3810–3829,, 2012. a

Wright, P. J., Harper, J. T., Humphrey, N. F., and Meierbachtol, T. W.: Measured basal water pressure variability of the western Greenland Ice Sheet: Implications for hydraulic potential, J. Geophys. Res.-Earth Surf., 121, 1134–1147,, 2016. a

Short summary
During glacial times in Switzerland, glaciers of the Alps excavated valleys in low-lying regions that were later filled with sediment or water. How glaciers eroded these valleys is not well understood because erosion occurred near ice margins where ice moved slowly and was present for short times. Erosion is linked to the speed of ice and to water flowing under it. Here we present a model that estimates the location of water channels beneath the ice and links these locations to zones of erosion.