Journal cover Journal topic
The Cryosphere An interactive open-access journal of the European Geosciences Union
Journal topic
The Cryosphere, 12, 2515–2544, 2018
https://doi.org/10.5194/tc-12-2515-2018
The Cryosphere, 12, 2515–2544, 2018
https://doi.org/10.5194/tc-12-2515-2018

Research article 07 Aug 2018

Research article | 07 Aug 2018

# Numerical reconstructions of the flow and basal conditions of the Rhine glacier, European Central Alps, at the Last Glacial Maximum

Numerical reconstructions of the flow and basal conditions of the Rhine glacier, European Central Alps, at the Last Glacial Maximum
Denis Cohen1, Fabien Gillet-Chaulet2, Wilfried Haeberli3, Horst Machguth3,4, and Urs H. Fischer5 Denis Cohen et al.
• 1Department of Earth and Environmental Science, New Mexico Tech, Socorro, NM, USA
• 2Univ. Grenoble Alpes, CNRS, IRD, Grenoble INP, IGE, Grenoble, France
• 3Department of Geography, University of Zurich, Zurich, Switzerland
• 4Department of Geosciences, University of Fribourg, Fribourg, Switzerland
• 5Nagra, Wettingen, Switzerland

Correspondence: Denis Cohen (denis.cohen@gmail.com)

Abstract

At the Last Glacial Maximum (LGM), the Rhine glacier in the Swiss Alps covered an area of about 16 000 km2. As part of an integrative study about the safety of repositories for radioactive waste under ice age conditions in Switzerland, we modeled the Rhine glacier using a thermodynamically coupled three-dimensional, transient Stokes flow and heat transport model down to a horizontal resolution of about 500 m. The accumulation and ablation gradients that roughly reproduced the geomorphic reconstructions of glacial extent and ice thickness suggested extremely cold (${T}_{\text{July}}\sim \mathrm{0}\phantom{\rule{0.125em}{0ex}}{}^{\circ }\mathrm{C}$ at the glacier terminus) and dry (∼10 % to 20 % of today's precipitation) climatic conditions. Forcing the numerical simulations with warmer and wetter conditions that better matched LGM climate proxy records yielded a glacier on average 500 to 700 m thicker than geomorphic reconstructions. Mass balance gradients also controlled ice velocities, fluxes, and sliding speeds. These gradients, however, had only a small effect on basal conditions. All simulations indicated that basal ice reached the pressure melting point over much of the Rhine and Linth piedmont lobes, and also in the glacial valleys that fed these lobes. Only the outer margin of the lobes, bedrock highs beneath the lobes, and Alpine valleys at high elevations in the accumulation zone remained cold based. The Rhine glacier was thus polythermal. Sliding speed estimated with a linear sliding rule ranged from 20 to 100 m a−1 in the lobes and 50 to 250 m a−1 in Alpine valleys. Velocity ratios (sliding to surface speeds) were >80 % in lobes and ∼60 % in valleys. Basal shear stress was very low in the lobes (0.03–0.1 MPa) and much higher in Alpine valleys (>0.2MPa). In these valleys, viscous strain heating was a dominant source of heat, particularly when shear rates in the ice increased due to flow constrictions, confluences, or flow past large bedrock obstacles, contributing locally up to several watts per square meter but on average 0.03 to 0.2 W m−2. Basal friction acted as a heat source at the bed of about 0.02 W m−2, 4 to 6 times less than the geothermal heat flow which is locally high (up to 0.12 W m−2). In the lobes, despite low surface slopes and low basal shear stresses, sliding dictated main fluxes of ice, which closely followed bedrock topography: ice was channeled in between bedrock highs along troughs, some of which coincided with glacially eroded overdeepenings. These sliding conditions may have favored glacial erosion by abrasion and quarrying. Our results confirmed general earlier findings but provided more insights into the detailed flow and basal conditions of the Rhine glacier at the LGM. Our model results suggested that the trimline could have been buried by a significant thickness of cold ice. These findings have significant implications for interpreting trimlines in the Alps and for our understanding of ice–climate interactions.

1 Introduction

Figure 1Map of Swiss Alps at the LGM showing maximum extent of ice cover (from Bini et al.2009) with outline of Rhine glacier basin in blue. Source: Bundesamt für Landestopografie swisstopo. Outlines of Figs. 2 and 3 are shown in red.

During the Last Glacial Maximum (LGM), the Alps were heavily glaciated and the Rhine glacier formed a large transection glacial complex that drained ice from the Alps north into the central Swiss and southwestern German Alpine forelands (Fig. 1). Repeated glacial advances into the lowlands during the Pleistocene sculpted the present-day landscape, forming emblematic valleys, horns, and arêtes in the Alps; deep, narrow lakes covering glacially excavated overdeepenings partially filled with glacial deposits; and moraines, outwash planes, terraces, and other depositional landforms in the lowlands. Numerous geomorphic studies since the beginning of the twentieth century have helped constrain the geometry and some flow characteristics of the Rhine glacier at its maximum extent during the last glaciation. Geomorphic mappings of terminal and lateral moraines delineated the extent of ice advances at glacial maxima and at various intermediate positions . Erratics, tills, and other sediment deposits provided information about flow paths and provenances . Presence of till indicated depositional environments while evidences of quarried bedrock, glacial polish, and striated rock surfaces identified areas with predominantly temperate (at the pressure melting point of ice), wet-based, erosive basal conditions . In the accumulation area, trimlines represent the upper limit of glacial erosion and define the minimum elevation of the ice surface . This large amount of geomorphic information for the Rhine glacier, probably the best documented of any paleo-glacier, has been used to create detailed maps of reconstructed ice extent and ice surface elevation at the LGM (e.g., ; ; ; ; ).

These detailed maps have been used to infer quantitative glaciological characteristics of the Rhine glacier at the LGM. In the ablation area, the first paleo-glaciological studies using simplified two-dimensional models indicated that thin, flat, and extended lobes of the large piedmont glaciers spreading out over much of the Swiss Plateau were polythermal with temperate basal conditions and low driving stresses (about 0.03 MPa, Haeberli and Penz1985). At the LGM, central Europe experienced extremely cold and dry conditions with the penetration of winter sea ice to low latitudes in the Atlantic Ocean. Correspondingly, the closure of the primary humidity source north of the Alps implied that most of the moisture feeding glaciers in the Alps had a southern, Mediterranean origin . Numerical models of LGM paleo-climate indicate a dry continental northern and northwestern Europe with more meridional storm tracks in comparison to present-day conditions. As a result, glaciers north of the ice divide had small ice throughflow with small surface velocities. In the northern lobes of piedmont glaciers, estimated average surface velocity was only about 25 m a−1 for the Rhine lobe . Average ice flow velocity through the main outlets of the Rhine glacier was estimated to be less than 200 m a−1 , a value relatively small compared to present-day ice bodies of comparable size. Extreme cold conditions also favored the development of permafrost up to about 150 m thick north of the marginal zone of the Rhine glacier . Subsurface temperatures and groundwater flow conditions must have been strongly influenced by the presence of extended surface and subsurface ice (Speck1994). Under these conditions, the glacier margins likely consisted of cold ice frozen to the subglacial permafrost, with limited basal sliding and basal meltwater flow, conditions that would have prevented significant glacial erosion there. Yet, erosional features such as overdeepenings are found close to earlier positions of the ice margin of the Rhine glacier (Fig. 2Preusser et al.2011; Dehnert et al.2012), necessitating basal conditions favorable for erosion. These could have occurred as a result of more humid conditions with higher ice flow velocity and increased sliding during ice advances across the Swiss Plateau, or during the rapid down-wasting of the ice mass that is likely to have taken place during the retreat phases, generating large quantities of water necessary for subglacial erosion. Another possibility is excavation during earlier ice ages with larger ice extents and warm-based conditions at sites with cold and frozen ice margins during the LGM. Alternatively, these early simple models may not have captured the full complexity of basal conditions near the glacier terminus at the LGM. Early two-dimensional models of the Rhine glacier at the LGM were not able to include the complexity of ice flow inherent to transection glaciers. Simple one-dimensional gravity-driven flow approximations of flow velocities from three-dimensional ice surface reconstructions (e.g., Benz-Meier2003) provided interesting first-order results but were physically not consistent when sliding was an important component of flow, neglected important stresses (longitudinal and transverse) relevant to the complex ice flow patterns of transection glaciers , and ignored glacier thermodynamics. An estimate of potential erosion for the Rhine glacier using this approximation did not find a correlation between ice sliding speed and areas of overdeepenings in the Rhine lobe.

Because of uncertainties in the reconstructed ice surface geometry and derived glaciological quantities in both the accumulation zone (more limited field evidences, trimline uncertainties, poorly known accumulation at the LGM) and in the ablation zone (poorly constrained temperate versus frozen basal conditions, ice flux, basal shear stress, and melt rates), geomorphic reconstructions of ice surface geometry should be verified against a three-dimensional ice flow model of the Rhine glacier. Flow conditions (patterns, velocities, stresses) obtained from geomorphic ice surface reconstruction should abide by the fundamental principles of glacier dynamics . For example surface flow directions should follow major ice drainages and computed vertical surface velocity from reconstructed ice surface should match mass balance conditions imposed by the LGM climate (e.g., Haeberli1991). These conditions may not always be satisfied with geomorphically reconstructed ice surfaces and warrant alternative approaches for paleo-reconstructions of ice bodies. The limits of geomorphic reconstructions and of simple stress-driven calculations necessitate a numerical modeling approach that can more precisely estimate the flow of the Rhine glacier at the LGM.

Finally, knowing the areas of temperate basal conditions, where sliding dominates, is important for addressing the safety of deep geological repositories for radioactive wastes. The long-term management of these wastes produced through the use of radioactive materials in power production, industry, research, and medicine entails their containment and isolation for over hundreds of thousands of years. Over such extended time periods, the performance of repositories in midlatitude and high-latitude regions can be affected by impacts from future ice age conditions. Several countries have developed programs to investigate potential future ice-related environmental changes and their effects at depth . The main concerns are deep erosion by glaciers or ice sheets, penetration of permafrost to great depth, and changes in groundwater hydrology due to permafrost and ice loading and their complex interactions. To address these issues, two main sources of information are being used: (i) qualitative and quantitative information about regional climate and ice conditions during past ice ages, and especially the LGM, as interpreted from paleo-climatic and paleo-ice proxies, and (ii) numerical modeling of complex and strongly coupled ice–climate systems. The present study was carried out within the framework of considerations concerning the long-term safety of repositories for radioactive waste in northern Switzerland (see Fig. 3).

Here, we use a fully three-dimensional, numerical, thermomechanical ice flow model that solves Stokes equations to investigate in detail the general characteristics of the Rhine glacier and to critically reflect on the accuracy of the geomorphic reconstructions with respect to ice flow physics and LGM climate. The computational burden of solving Stokes equations implies that short (a few thousand years at most) transient simulations around the LGM are used only to seek steady-state solutions, neglecting transient effects of climate and other processes (e.g., isostatic adjustment). Stokes equations are solved using Elmer/Ice , an open-source finite-element code for ice flow. The ice flow model is driven by a simple mass balance model parameterized by two mass balance gradients, one for the accumulation zone, one for the ablation zone, and a given equilibrium line altitude (ELA). Parameterization of temperature is based on a given temperature at the ELA and a lapse rate. The model yields the full three-dimensional velocity and temperature fields; details of surface, englacial, and sliding speeds; basal temperatures; and shear stresses from which ice fluxes, flow patterns, and areas of temperate basal condition can be derived.

Figure 2Depth of overdeepened valleys in the Rhine glacier basin . Outlines of maximum glacial extent at the LGM are shown in blue.

2 The Rhine glacier model

We model the flow of the Rhine glacier using the full equations of mass and momentum, the Stokes flow equations, coupled to the heat equation over the Rhine glacier basin (Fig. 1) at a horizontal resolution of about 500 m using the finite-element, open-source code Elmer/Ice . Stokes flow is more accurate than shallow ice and shallow shelf approximation in zones in which gradients in stresses are significant, like in complex transection glaciers found in Alpine settings and where sliding is also significant . Because of computational cost, transient simulations are limited to several thousand years at most, seeking steady-state solutions for a constant LGM climate (temperature, mass balance gradients). Our domain of computation for the Rhine glacier includes all basins that drain into present-day Lake Constance and Lake Zurich. These basins straddle four countries: Switzerland, Germany, Liechtenstein, and Austria. Present-day topographic divides with the Rhone basin to the west, the Ticino basin to the south, and the Inn basin to the southeast delineate sources of ice in the accumulation area. For the ablation area to the north, the maximum glacial extent at the LGM is extended about 50 km northward of the Rhine and Linth lobes to delimit the extent of the model. This allows the Rhine glacier to advance further north than its LGM extent in the numerical simulations if necessary. On the western side, the divide between the Reuss and Linth lobes serves as the model boundary. The eastern limit of the Rhine lobe is used to delineate the eastern edge of the model.

The basal topography of our ice flow model uses the present-day topography (initially obtained from a 25 m resolution digital elevation model). This topography was modified to remove present-day ice thickness using published results of an inversion model by Linsbauer and colleagues and lake bathymetry where available (Lake Constance and Lake Zurich; Benz-Meier2003) and by depressing the topography to reproduce isostatic adjustments at the LGM (up to 130 m) using the model of . Figure 3 shows the basal map used for the model.

Figure 3Basal topography (from Benz-Meier2003) used in all simulations with names of major valleys, present-day lakes, cities, and mountains mentioned in text. The blue outline of the Rhine glacier at the LGM is from . See text for details. The three proposed siting regions for high-level waste repositories in Switzerland are shown in red. All three sites are located within the ice extent of earlier glaciations larger than the LGM.

For the ice surface, we use either the digitally available reconstructed glacier surface of or an ice surface obtained from a previous numerical simulation. The reconstruction of is based on the earlier works of and also of . We refer to it as despite its earlier origin. A two-dimensional, fixed, unstructured triangular mesh over the domain is created using Gmsh and then extruded into hexahedral elements in the vertical direction between the basal surface and the ice surface.

## 2.1 Thermomechanical model

Stokes equations describe the mass and momentum balance for a viscous fluid. Together with the equation for energy, this system forms a thermomechanical problem of five coupled scalar equations with five field variables to solve for p the pressure, v the velocity vector (u, v, w), and T the temperature. In addition, for ice flow, the glacier surface elevation, ζ, is unknown and requires an additional equation, the kinematic boundary condition.

## 2.2 Stokes equations

The conservation of mass for an incompressible fluid (divergence-free velocity field) is

$\begin{array}{}\text{(1)}& \mathrm{\nabla }\cdot \mathbit{v}=\mathrm{0},\end{array}$

where ∇⋅ indicates the divergence. The conservation of momentum is

$\begin{array}{}\text{(2)}& -\mathrm{\nabla }p+\mathit{\eta }{\mathrm{\nabla }}^{\mathrm{2}}\mathbit{v}+\left[\mathrm{\nabla }\mathbit{v}+{\left(\mathrm{\nabla }\mathbit{v}\right)}^{T}\right]\mathrm{\nabla }\mathit{\eta }+\mathit{\rho }\mathbit{g}=\mathrm{0},\end{array}$

where g is the vector of gravity, ρ is the density of ice, indicates the gradient, and η is the ice viscosity given by Glen's flow law

$\begin{array}{}\text{(3)}& \mathit{\eta }=\left\{\begin{array}{l}\frac{\mathrm{1}}{\mathrm{2}}{\left(EA\left(T\right)\right)}^{-\mathrm{1}/n}{d}^{-\left(\mathrm{1}-\mathrm{1}/n\right)},d>{d}_{o},\\ \\ \frac{\mathrm{1}}{\mathrm{2}}{\left(EA\left(T\right)\right)}^{-\mathrm{1}/n}{d}_{o}^{-\left(\mathrm{1}-\mathrm{1}/n\right)},d\le {d}_{o},\end{array}\right\\end{array}$

where $d=\sqrt{\left(\frac{\mathrm{1}}{\mathrm{2}}\text{tr}\left({D}^{\mathrm{2}}\right)\right)}$ is the effective strain rate, do is the critical strain rate, $D=\frac{\mathrm{1}}{\mathrm{2}}\left[\left(\mathrm{\nabla }\mathbit{v}\right)+{\left(\mathrm{\nabla }\mathbit{v}\right)}^{T}\right]$ the strain rate tensor, n the power-law exponent, T the temperature, A(T) the rate factor, and E the flow-enhancement factor taken equal to 1 in all simulations. An Arrhenius relationship is used for A,

$\begin{array}{}\text{(4)}& A\left(T\right)={A}_{o}\left(T\right)\mathrm{exp}\left(\frac{-Q\left(T\right)}{R\phantom{\rule{0.125em}{0ex}}T}\right),\end{array}$

where Ao is the pre-factor, Q is the activation energy, and R is the gas constant. Different values of Ao and Q are used depending on whether T is greater or less than −10C.

## 2.3 The heat equation

The heat equation for a viscous fluid is

$\begin{array}{}\text{(5)}& \mathit{\rho }c\left(T\right)\left(\frac{\partial T}{\partial t}+\mathbit{v}\phantom{\rule{0.125em}{0ex}}\mathrm{\nabla }T\right)=\mathrm{\nabla }\left(\mathit{\kappa }\left(T\right)\mathrm{\nabla }T\right)+\mathrm{4}\mathit{\eta }{d}^{\mathrm{2}},\end{array}$

where c(T) is the heat capacity of ice and κ(T) the heat diffusivity of ice, both functions of temperature. Parameterization of these equations and values of model parameters for the ice flow model are given in Table 1. The last term in the equation above describes viscous dissipation.

Table 1Model parameter values.

## 2.4 Basal boundary conditions

### 2.4.1 Geothermal heat flow

Geothermal heat flow in the area of the LGM Rhine glacier is highly variable, ranging from 0.06 to 0.12 W m−2 with low values in the high Alps and high values in geothermally active regions of the Swiss Plateau. We use the heat flow data of obtained from numerous measurements of temperature gradients in boreholes and of rock thermal properties. Present observed temperature gradients, however, are affected by topographic effects, long-term changes in climate and erosion, groundwater flow, and past ice sheet cover in a complex, nonlinear way. Correcting for these effects is difficult and is not performed in ; neither is it carried out here. It would require, at a minimum, imposing the heat flux at the base of a 2 to 3 km thick bedrock below the ice and hence would necessitate solving the heat equation in the bedrock beneath the ice, a task beyond the scope of the present study.

### 2.4.2 Basal sliding

At the bed, ice can slide over bedrock when water is present, which occurs when the temperature at the bed equals the melting temperature (so-called temperate bed). Subfreezing sliding can also occur (e.g., Shreve1984) but is small and negligible for ice fluxes or erosion for the timescale considered here (hundreds to thousands of years). It is neglected here. A commonly used approach for the sliding speed is to link it to the basal shear stress, i.e.,

$\begin{array}{}\text{(6)}& {\mathbit{\tau }}_{\mathrm{b}}=C\phantom{\rule{0.125em}{0ex}}‖{\mathbit{v}}_{\mathrm{s}}{‖}^{m-\mathrm{1}}{\mathbit{v}}_{\mathrm{s}},\end{array}$

where vs is the sliding speed vector, τb the basal shear stress vector, m an exponent, and C a constant that encapsulates the effects of water pressure, bed roughness, etc. The parameter m usually takes values between 1 and 1∕3. In all simulations we assume m=1 (linear sliding), acknowledging that this does not realistically model several modes of basal ice motion such as flow over till layers with subglacial water flow couplings (e.g., weak ice–bed couplings; Iverson2010), flow over bedrock protuberances with ice–bed separation and water-filled cavities , or fast ice sheet flows (e.g., Engelhardt and Kamb1998). Given the large uncertainties in predicting and modeling glacial sliding , our approach for this paleo-ice flow model is to keep the sliding law simple.

Ice sliding over the substrate occurs when temperatures reach the melting temperature. Below this temperature, vs=0. To simulate this behavior we make the sliding coefficient C in the sliding rule above a function of temperature T ,

$\begin{array}{}\text{(7)}& C\left(T\right)=\left({C}_{o}-{C}_{\mathrm{1}}\right)\mathrm{exp}\left[-\frac{{T}^{\prime }}{\mathit{\gamma }}\right]+{C}_{\mathrm{1}},\end{array}$

where Co and C1 are the sliding coefficients for temperate and cold conditions, respectively, and ${T}^{\prime }={T}_{\text{pmp}}-T$, where Tpmp is the temperature at the pressure melting point. The parameter γ, sometimes called the sub-melt sliding parameter , adjusts the range of temperature over which this transition occurs. For numerical stability we use γ=2.

### 2.4.3 Thermal boundary condition

The boundary condition for the basal temperature is

$\begin{array}{}\text{(8)}& \mathit{\kappa }\left(T\right)\mathrm{\nabla }T\cdot \mathbit{n}={q}_{\text{geo}}-\mathbit{t}\cdot {\mathbit{v}}_{\mathrm{s}},\end{array}$

where n is the normal to the glacier bed, qgeo is the geothermal heat flux, and t is the stress vector at the glacier bed. The last term of the equation is the thermal energy due to basal friction. In most simulations this term is assumed to be zero. When the glacier bed is at the pressure melting point, Eq. (8) no longer holds but the difference between the left and right sides of this equation can be used to compute the melt rate.

## 2.5 Boundary conditions at the ice surface

Solving for ice velocity, temperature, and the elevation of the ice surface requires three types of boundary conditions at the ice surface: (i) a kinematic boundary condition that describes the motion of the ice surface as a function of accumulation (or ablation) and ice flow, (ii) an expression of the stress-free condition at the ice surface, and (iii) the ice temperature.

### 2.5.1 Stress-free surface

The glacier surface is a free surface in contact with the atmosphere. This surface supports no shear stress. This is simply expressed as $\mathbit{t}\cdot \mathbit{n}=\mathrm{0}$, where t is the stress vector at the ice surface and n is the unit normal pointing outward.

### 2.5.2 Kinematic boundary condition

The ice surface moves up or down depending on the balance between mass flux across the glacier surface and the divergence of the velocity field. This is expressed mathematically by the kinematic boundary condition

$\begin{array}{}\text{(9)}& \frac{\partial \mathit{\zeta }}{\partial t}=\stackrel{\mathrm{˙}}{b}+w-u\frac{\partial \mathit{\zeta }}{\partial x}-v\frac{\partial \mathit{\zeta }}{\partial y},\end{array}$

where ζ is the glacier surface elevation and $\stackrel{\mathrm{˙}}{b}$ is the specific balance rate (specific balance for short), i.e., the volumetric mass flux of ice per unit time across the glacier surface, the accumulation and ablation function.

### 2.5.3 Surface mass balance

For the accumulation–ablation function ($\stackrel{\mathrm{˙}}{b}$ in Eq. 9), we choose a simple parameterization represented by two mass balance gradients, one for the accumulation area and one for the ablation area, and a maximum threshold value for the maximum accumulation,

$\begin{array}{}\text{(10)}& \stackrel{\mathrm{˙}}{b}=\left\{\begin{array}{lr}min\left({\stackrel{\mathrm{˙}}{b}}_{\text{acc}}^{\text{max}},{\mathit{\beta }}_{\text{acc}}\left(z-{z}_{\text{ela}}\right)\right),& z\ge {z}_{\text{ela}},\\ {\mathit{\beta }}_{\text{abl}}\left(z-{z}_{\text{ela}}\right),& z<{z}_{\text{ela}},\end{array}\right\\end{array}$

where z is the elevation of the ice surface , βacc and βabl are the accumulation and ablation mass balance gradients, respectively, and ${\stackrel{\mathrm{˙}}{b}}_{\text{acc}}^{\text{max}}$ is an upper bound for the accumulation rate.

This parameterization is a simplification of the actual and mostly unknown spatial patterns of accumulation and ablation processes at the LGM. Together with an equation for the surface temperature (Sect. 2.5.4), we have effectively decoupled the mass balance from the energy. A rate of ablation based on the number of positive degree days (PDD model) would have been more physical. However, because of large uncertainties in LGM climate (including annual temperature amplitudes), previous applications of the PDD approach had to rely on present-day temperature distribution minus an offset . Furthermore, PDD factors needed to compute surface melt rates are known to vary substantially and choosing suitable factors for LGM conditions is challenging. The same is true for accumulation, which requires knowing patterns of precipitation and temperature. Since LGM climate is known to have behaved rather differently from today owing to the southward displacement of the atmospheric polar front in the North Atlantic , it is questionable whether the added complexity of a PDD model results in a more accurate representation of LGM mass balance distribution.

Table 2 lists the chosen mass balance gradients. Different values of βacc and ${\stackrel{\mathrm{˙}}{b}}_{\text{max}}$ were selected to represent a range of dry to wet climates while values of βabl were chosen to be in agreement with either earlier estimates for the LGM Rhine glacier or present-day Greenland values measured in regions where temperatures and precipitation are similar to estimated LGM conditions in northern Switzerland.

### 2.5.4 Surface temperature

Because of the uncoupling of mass balance and energy, surface temperature only influences ice temperature and rheology (see Eq. 3) but not the rate of accumulation or ablation. As carried out in many ice modeling studies , we assume that the ice surface temperature is equal to the mean annual air temperature. Based on earlier climate reconstructions we assume that the mean annual temperature at the equilibrium line is −12C, a value within a range of estimations (−15 to −10C; Haeberli1983, 2010). We assume that changes in surface temperature with elevation are linearly related to a lapse rate, γa,

$\begin{array}{}\text{(11)}& {T}_{\text{surf}}\left(z\right)={T}_{\text{ela}}+{\mathit{\gamma }}_{a}\left(z-{z}_{\text{ela}}\right),\end{array}$

where Tela and zela are the temperature and elevation at the equilibrium line, respectively. In all model simulations ${\mathit{\gamma }}_{a}=-\mathrm{6}$${}^{\circ }\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{km}}^{-\mathrm{1}}$, an intermediate value based on an estimate by (−7${}^{\circ }\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{km}}^{-\mathrm{1}}$) for the Rhine glacier at the LGM and contemporary lapse rates for polar regions (−5 to −4${}^{\circ }\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{km}}^{-\mathrm{1}}$Marshall et al.2007; Machguth and Cohenunpublished). Estimates of the equilibrium line altitude at the LGM vary from 944 m to 1200 m . Benz's estimate is based on the assumption that the accumulation area ratio (the ratio of the accumulation area of the glacier to the total surface area) is 2∕3. Another estimate by using the same method but a slightly different ice surface reconstruction and hypsographic curve found an ELA of 1000 m. In our simulation, ELA ranges from 1000 to 1200 m (see Table 2).

## 2.6 Lateral boundaries

Nodes on lateral boundaries are either no-flow nodes, for which the normal component of velocity perpendicular to the boundary is zero, or outflow nodes, for which the natural stress-free boundary condition is applied. Outflow nodes are used on the northern boundary. No-flow nodes are at the glacier basin boundary in the accumulation zone. The horizontal temperature gradient is assumed to be zero (no-flux condition) along all lateral boundaries.

## 2.7 Temperature initialization

One of the most difficult fields to initialize is the englacial temperature (velocities and pressure are easily computed from the Stokes equations once temperature is known). In reality, the temperature of an ice sheet depends on the interplay among climate, flow, and geothermal heating so that the temperature at any instant should require, in theory, knowledge of the full flow and climate history that affected the ice sheet. For the problem at hand, this is not known. We initialize the temperature field in the ice assuming the temperature depends on an accumulation or ablation rate (vertical advection) and is independent of horizontal advection. The advantage of these assumptions is that an analytical expression for temperature exists from (see Cuffey and Paterson2010, p. 410). The disadvantage of this initialization is that horizontal advection is neglected and thus the temperature distribution is out of equilibrium and simulations necessitate a long spin-up time.

The vertical temperature distribution for an ice thickness between 0 and H with vertical advection given by $w=-\stackrel{\mathrm{˙}}{b}\phantom{\rule{0.125em}{0ex}}z/H$ is, for the accumulation zone ($\stackrel{\mathrm{˙}}{b}>\mathrm{0}$),

$\begin{array}{}\text{(12)}& T={T}_{\text{surf}}+{z}_{\star }\frac{\sqrt{\mathit{\pi }}}{\mathrm{2}}{\left[\frac{\mathrm{d}T}{\mathrm{d}z}\right]}_{\text{bed}}\left[\text{erf}\left(z/{z}_{\star }\right)-\text{erf}\left(H/{z}_{\star }\right)\right],\end{array}$

where [dT∕dz]bed is the gradient of temperature at the bed and ${z}_{\star }^{\mathrm{2}}=\mathrm{2}{\mathit{\alpha }}_{T}H/\stackrel{\mathrm{˙}}{b}$. For ablation ($\stackrel{\mathrm{˙}}{b}<\mathrm{0}$), the temperature distribution is

$\begin{array}{ll}T={T}_{\text{surf}}& +{\left[\frac{\mathrm{d}T}{\mathrm{d}z}\right]}_{\text{bed}}\left\{{z}_{\star }\mathrm{exp}\left({z}^{\mathrm{2}}/{z}_{\star }^{\mathrm{2}}\right)\phantom{\rule{0.125em}{0ex}}\mathcal{D}\left(z/{z}_{\star }\right)\right\\\ \text{(13)}& & -{z}_{\star }\mathrm{exp}\left({H}^{\mathrm{2}}/{z}_{\star }^{\mathrm{2}}\right)\phantom{\rule{0.125em}{0ex}}\mathcal{D}\left(H/{z}_{\star }\right)},\end{array}$

where 𝒟() is the Dawson integral

$\begin{array}{}\text{(14)}& \mathcal{D}\left(x\right)={e}^{-{x}^{\mathrm{2}}}\underset{\mathrm{0}}{\overset{x}{\int }}{e}^{{t}^{\mathrm{2}}}\mathrm{d}t.\end{array}$

Using the accumulation–ablation function given in Eq. (10), the ice surface temperature, and a given heat flux, the temperature distribution in the accumulation and ablation zones is computed everywhere in the glacier using Eqs. (12) and (13). The computed temperature sometimes exceeds the melting temperature when the ice is thick. This is not surprising since nowhere in the model of is there information about a maximum temperature or about a phase transition. Limiting the temperature to the melting temperature is performed as a post-processing step. This obviously breaks apart the conservation of energy but is a quick fix to obtain a first-order approximation of temperature. Temperatures warmer than the melting temperature are obtained because, in the accumulation zone, thick ice insulates ice at depth from the cold atmospheric conditions and geothermal heat flux adds heat to the ice from the bottom. In the ablation zone, the process is similar, but in addition, temperate ice is advected upward.

Table 2Summary of numerical (this study) and selected geomorphic reconstructions. zela is the elevation of the equilibrium line. βabl and βacc are the ablation and accumulation mass balance gradients, respectively. ${\stackrel{\mathrm{˙}}{b}}_{\text{acc}}^{\text{max}}$ is the maximum accumulation rate. ${\stackrel{\mathrm{‾}}{b}}_{\text{acc}}$ is the average accumulation rate. ${\stackrel{\mathrm{‾}}{b}}_{\text{net}}$ is the glacier net mass balance. ${\stackrel{\mathrm{˙}}{b}}_{\text{term}}$ is the ablation rate at the terminus. AAR is the accumulation area ratio. A is the glacier area and V its volume. ts is the simulated time. Numerical simulations s1 through s5 are ranked from driest to wettest based on the value of the average accumulation rate, ${\stackrel{\mathrm{‾}}{b}}_{\text{acc}}$. Initial conditions for the ice surface are the geomorphic reconstruction of for s1, a simulation with βacc=0.1$\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{m}{\right)}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$, βabl=0.2$\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{m}{\right)}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$, and zELA=1200m that ran for 440 years using the reconstruction as the initial condition for s2, and a simulation with βacc=0.05$\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{m}{\right)}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$, βabl=0.2$\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{m}{\right)}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$, and zELA=900m that ran for 907 years using as the initial condition for s3, s4, and s5.

a. b. cEquilibrium is assumed.

## 2.8 Numerical simulations

Five simulations (Table 2), named s1 through s5, with different values of mass balance gradients (accumulation and ablation) and ELA, represent a range of climate scenarios from driest (s1) to wettest (s5). The range of ELAs correspond to the range estimated from previous reconstructions while mass balance gradients were chosen to loosely reproduce LGM climatic conditions. Because of computational cost, not all simulations could be run for over 3000 years due to differences in the rate of convergence of the nonlinear system of equations. Also, long computation time made a parametric study unfeasible. Instead, our choice of parameter values for these five simulations is based on sampling the parameter space using realistic values. Scenario s1 has the smallest accumulation rate and the smallest melt rate at the glacier terminus, representative of an extreme climate, perhaps colder and drier than LGM conditions. Simulation s5, which has the highest cutoff value for the maximum accumulation rate (Eq. 10), is meant to represent the other extreme: a wetter climate in the south that may have occurred as a result of moisture transport from south of the Alps that yielded significant accumulation in the high Alpine peaks but that remained relatively cold and dry in the north near the glacier terminus. Although our mass balance model does not include a north–south gradient, the increase in elevation southward combined with a higher cut-off value could conceptually represent such a directional gradient in mass balance. A possible alternative to creating a north–south gradient in mass balance would have been to use an ELA that decreased southward towards the source of moisture. Such a model was used by to fit observations of mass balance distribution for the Vatnajökull ice cap in Iceland. Here and in the remainder of this paper, the notion of colder–dryer versus warmer–wetter conditions is based solely on the values of mass balance in the accumulation and ablation zones, not on the temperature prescribed at the ice surface. Since in our model the ice surface temperature is decoupled from the mass balance, temperature has only lower-order effects on the ice flow dynamics, which is mainly controlled by the mass balance. Other parameters (basal topography, glaciological parameters) are identical for all simulations. All input parameters and key computed quantities for these simulations are summarized in Table 2. Note that not all simulations started with identical initial conditions. Since our goal is not a systematic study of the effect of parameter values on the Rhine glacier ice flow dynamics, we believe these changes in initial condition have little impact on the general results regarding ice flow conditions at the LGM. For comparison, the table also includes data for two geomorphic reconstructions: and .

Input and output quantities are different for the numerical and the geomorphic reconstructions. For the numerical simulations, model input parameters are the ELA (zela), the mass balance gradients (βabl and βacc), and the maximum rate of accumulation (${\stackrel{\mathrm{˙}}{b}}_{\text{max}}$). Model outputs are the average net accumulation (${\stackrel{\mathrm{‾}}{b}}_{\text{acc}}$), the average net glacier balance (${\stackrel{\mathrm{‾}}{b}}_{\text{net}}$), the specific balance rate at the terminus (${\stackrel{\mathrm{˙}}{b}}_{\text{term}}$) estimated at an elevation of 500 m, the accumulation area ratio (AAR, ratio of accumulation area to total area), and the glacier area (A) and its volume (V). Also indicated in the table is the simulated time, ts, ranging between 1500 and over 3000 years. For the geomorphic reconstructions, inputs and outputs are almost reversed. Inputs are the glacier area and volume obtained from field mapping and inferences of ice surface elevation contours, AAR, ${\stackrel{\mathrm{˙}}{b}}_{\text{term}}$, and ${\overline{b}}_{\text{net}}$ assumed to be zero (the glacier is in a state of dynamic equilibrium).

For the geomorphic reconstructions, the term ${\stackrel{\mathrm{˙}}{b}}_{\text{term}}$ is calculated from values of summer and mean annual LGM temperatures near the 500 m elevation level indicated in both and . These temperatures are converted to melt rates, first by estimating the number of PDDs using the model of Reeh (1991) and then by multiplying the PDD by a PDD factor, here taken equal to 6 mm PDD−1 , a value used for melting ice in Greenland. Higher PDD values (6–8 mm PDD−1) have been measured on mountain glacier (e.g., Braithwaite1995; Hock2003; Jóhannesson et al.1995) but LGM climate conditions in northern Switzerland may have likely been more similar to present-day Arctic conditions. Both and assume an AAR of 0.67 based on earlier studies of modern Alpine glaciers (e.g., Gross et al.1977) that assume zero net balance. This assumption, together with the glacier hypsometry, determines the ELA. Using the ELA and the melt rate at the terminus, one can compute the average ablation gradient, βabl. Assuming the glacier was at equilibrium, net accumulation equals net ablation. Net ablation can be calculated from the mass balance gradient computed with the method just described and the glacier hypsometry (area distribution of elevation) below the ELA. Using the same procedure but in reverse, the mass balance gradient in the accumulation area, βacc, can be calculated from the net accumulation (equals the net ablation) and the glacier hypsometry above the ELA. Numbers for the ablation and accumulation mass balance gradients of the geomorphic reconstructions calculated with this method are given in Table 2. Note that the average accumulation rate calculated for (0.50 m a−1, Table 2) is higher than the value cited in (0.30 m a−1) obtained from numbers cited in that assumed LGM precipitation was equal to 20 % of today's value (1.5 m a−1).

## 2.9 Numerical solutions

The open-source software Elmer/Ice is used to solve the set of equations and their boundary conditions using the finite-element method. The general three-step procedure for initializing and running simulations is as follows.

• Solve the transient Stokes flow equations together with the free surface for 50 years with a constant temperature field given by the initialization. The energy equation is turned off.

• Solve the steady-state energy equation only using the velocity field calculated at the end of the 50 years.

• Solve all fields simultaneously in transient mode (flow, energy, free surface) for several thousand years, possibly reaching a steady state. Some simulations were stopped earlier due to computational time constraints.

3 Model results

## 3.1 Ice extent, surface elevation, and ice thickness

Figure 4Ice surface elevation in meters above sea level for (a)  and (b–f) simulations s1 through s5. Black contours are every 500 m. The 2500 m contour is shown in yellow. The magenta contour indicates the equilibrium line altitude. The equilibrium line is 1200 m except for simulations s2 (1000 m) and s4 (1100 m). Topography is shown in ice-free areas with a different color scheme. In (a) the northern extent of the Rhine and Linth glacier lobes are based on mapping of terminal moraines.

Ice extent and ice surface elevation of simulations s1 through s5 are shown in Fig. 4. In almost all subsequent figures that present numerical simulation results, the geomorphic reconstruction of or a numerical simulation based on his ice surface elevation is included for comparison (e.g., Fig. 4a).

Out of the five simulations shown in Fig. 4, simulation s1 appears to have essentially reached steady state (net glacier balance ${\stackrel{\mathrm{˙}}{b}}_{\text{net}}=\mathrm{0.02}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$; see Table 2) and simulation s3 seems close to it (${\stackrel{\mathrm{˙}}{b}}_{\text{net}}=\mathrm{0.19}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$). Whether these states are true equilibrium states is debatable as negative feedbacks associated with topography and changes in mass balance with elevation can lead to complete ice sheet collapse in simple models . Whether this applies to our model here is unknown. Other simulations remained out of equilibrium conditions during the entire simulation time with the ice mass still increasing. The two simulations that are closest to an apparent steady state (simulations s1 and s3) have an ice surface area smaller than the LGM extent defined by moraines but a margin shape and configuration that resembles the outline of the LGM terminal moraines (e.g., ice-free Hörnli ridge in between the Rhine and Linth lobes, dendritic outline of ice extent). Of these two simulations, simulation s3 is closer to the LGM margin with an ice surface area of 14 553 km2 while simulation s1 is only 12 389 km2 (the geomorphic reconstructions yielded values around 16 000 km2). The ice surface area of simulation s5 is similar to the geomorphic reconstruction (16 400 km2). Simulation s4 is slightly larger (17 706 km2); s2 is slightly smaller (15 813 km2). Their configurations, however, are significantly different, with most Alpine peaks and the Hörnli ridge covered in ice, and with an ice margin that is less dendritic. Note that these three simulations (s2, s4, and s5) have not yet reached steady state and have a net positive mass balance (Table 2) with an ice mass that is still growing and a margin that continues to move north, beyond the LGM margin.

Only the modeled ice surface elevation of simulation s1 resembles the geomorphic reconstruction of (Fig. 4a, b). The 2500 m elevation contours (shown in yellow in Fig. 4) of the geomorphic reconstruction and of simulation s1 in the accumulation area are relatively close (within a few kilometers). Also, the 1200 m contour lines (the ELA for simulation s1) have very similar positions across the piedmont glacier that forms the Linth lobe (the location where the glacier flows through its narrowest point before spreading onto the lowlands), while for the Rhine lobe the 1200 m contour line bulges out slightly into the foreland in the numerical simulation. The ice surface elevations of simulations s2–s5 are significantly higher than the geomorphic reconstruction. In simulations s2–s4, the 3000 m elevation contour is found roughly at the same location as the 2500 m contour in the geomorphic reconstruction. Simulation s5 has a higher ice surface with the 3200 m elevation contour at roughly the same position as the 2500 m contour of the geomorphic reconstruction. This indicates an ice thickness 500 to 700 m higher for simulations s2–s5 than in the geomorphic reconstruction or in simulation s1. Ice thickness for all numerical simulations is shown in Fig. 5. Maximum ice thickness in the geomorphic reconstruction of is slightly greater than 1500 m in the Vorderrhein (see Fig. 3 for location). This same 1500 m contour has only slightly expanded in simulation s1. However, in all other simulations, maximum ice thickness exceeds 2000 m in the same area.

Figure 5Ice thickness in meters for (a)  and (b–f) simulations s1–s5.

## 3.2 Ice surface and hypsometric curve

Hypsometric curves are useful to show the area distribution of elevations of the ice surface. Figure 6 shows hypsometric curves for the five numerical simulations and for the and geomorphic reconstructions. The curves for the numerical simulations are shifted to the right in comparison to those for the geomorphic reconstructions, indicating more ice at higher altitudes. Only simulation s1 converges to the geomorphic reconstructions for high elevations (above 2500 m). For the numerical simulations, the values of the AAR can be read directly from the graph, the AAR being the intersection of the ELA with the hypsometric curve. All numerical reconstructions have a higher AAR than the value of 0.67 assumed for the two geomorphic reconstructions. For simulations s1 and s3, which are close to steady state, their AAR values are 0.69 and 0.71, respectively. Simulations s2, s4, and s5, which have not reached steady state, have higher AAR values: 0.88, 0.70, and 0.76, respectively (see Table 2). Since these simulations are still in a phase of ice growth, their steady-state AAR values are likely to be smaller. Note also that the lower the ELA, the higher the AAR.

Figure 6Glacier hypsometric curves for geomorphic and numerical reconstructions of the Rhine glacier system. Vertical colored lines show the value of the AAR on the x axis for each reconstruction.

## 3.3 Ice surface speed

Figure 7Ice surface speed in meters per year for (a)  (steady-state Stokes flow solution with no sliding and a uniform ice temperature of 0 C) and (b–d) simulations s1–s5. Also shown are ice flow lines separating major ice divides. Green: Linth–Walensee; black: Walensee–Rhine; orange: Rhine–Ill. See Sect. 4.3.2 for details and Fig. 1 for location of major glaciers. Note the different scale for speed between (a) and all other subfigures.

Figure 7 shows the computed ice surface velocities. All numerical simulations (s1 through s5, Fig. 7b–f) show patterns of glacial flow concentrated in the main trunks of Alpine valleys with the highest flow velocities in the Rhine glacier and at the confluence of the Linth and the Walensee glaciers. Flow velocities decrease rapidly as ice fans out into the Rhine lobe but radially oriented topographic lows that follow river drainage facilitate ice movement at higher speeds there. In contrast, calculated surface velocities for the geomorphically reconstructed ice surface show no such pattern of ice flow (Fig. 7a). Instead, flow is concentrated in zones of high slope gradients, giving rise to unrealistic patterns of ice surface speeds with no apparent channeling of ice through Alpine valleys. Also, ice is nearly stagnant in a large portion of the Rhine and Linth lobes, with surface velocities of less than 10 m a−1. For the geomorphically reconstructed simulation (Benz, Fig. 7a), the ice velocity was obtained from the solution of the Stokes flow equations in steady state with no sliding and with the reconstructed ice surface of fixed with the stress-free boundary condition. A uniform ice temperature of 0 C was also used to obtain the same rheological parameter (A in Eq. 3) as prescribed by . Although our Stokes flow model is different from the simple one-dimensional model of , our patterns of ice flow are comparable to his reconstruction (see Benz-Meier2003, Fig. 6.10) and ice velocities are comparable in magnitude for the most part (e.g., 2930 and 3130 m a−1 at the Sargans diffluence for our and Benz-Meier's reconstructions, respectively), although not for the areas with very high slope gradients near high Alpine peaks for which our model predicts velocities in excess of 5000 m a−1, far exceeding velocities predicted by Benz-Meier (∼2500m a−1). Table 3 summarizes maximum speeds for all simulations as well as other dynamical quantities described in the following sections.

Table 3Summary of ice dynamics quantities. Vmax is the maximum surface velocity. ${\stackrel{\mathrm{‾}}{V}}_{\text{lobe}}$ is the average surface velocity in the Rhine lobe. ${V}_{\mathrm{s}}^{\text{max}}$ is the maximum sliding velocity. ${\overline{\mathit{\tau }}}_{\text{bed}}$ is the average basal shear stress over the entire glacier bed. VH is the mean vertically integrated strain heating.

## 3.4 Basal conditions

Basal temperature and sliding speed are two key variables that describe basal conditions. Basal temperature relative to the melting point is shown in Fig. 8 as the difference between modeled ice temperature and the pressure melting point calculated as a function of ice pressure at the bed. The yellow contour outlines basal ice that is within 0.1 C of the pressure melting point, essentially temperate for all practical purposes. For the numerical simulation based on the geomorphic ice surface reconstruction (Fig. 8a), the bed is assumed temperate as in . For all other numerical reconstructions, large areas of the bed of Alpine glaciers and the majority of the Rhine and Linth lobes are at the melting temperature. Exceptions are the upper reaches of Alpine valleys and areas in the lobes where ice is thin and the cold temperatures at the ice surface diffuse down to the base of the glacier. Note also that, in general, the northeastern part of the Rhine lobe east of the Schussen valley (see Fig. 3 for location) is colder than the western part of the lobe, probably owing to rising hummocky topography in the east and to thinner ice cover there (see Fig. 5).

Figure 8Basal temperature relative to the melting point ($\mathrm{d}T=T-{T}_{\text{pmp}}$) in Kelvin for (a)  (temperate bed assumed) and (b–f) simulations s1–s5. Yellow contour indicates the extent of ice that is within 0.1 C of the pressure melting point.

Sliding speed is shown in Fig. 9. Following , zero basal velocity is assumed for the numerical geomorphic reconstruction (Fig. 9a). For other numerical simulations, the spatial patterns of sliding are directly related to the basal ice temperature (Eqs. 6 and 7). Significant sliding begins in the Rhine valley at the confluence of the Vorderrhein and the Hinterrhein with very limited sliding up-valley of this location. Downstream of the Vorderrhein–Hinterrhein confluence, temperate bed conditions prevail and sliding speeds reach their maximum values in the main trunk of the Rhine valley (Alpenrhein) around Sargans where the Rhine splits into the main Rhine and the Walensee glacier, and also in the lower portion of the Rhine valley at the Alpine gate before ice enters the German–Swiss forelands and the Lake Constance basin. Sliding is also significant at the confluence of the Linth and Walensee glaciers.

Figure 9Ice sliding speed in meters per year for (a)  (no sliding assumed) and (b–f) simulations s1–s5.

Following patterns of surface speed (Fig. 7), sliding speed is highest in the Alpenrhein in simulation s5 (Fig. 9f), reaching about 200 m a−1. Maximum sliding speed diminishes for simulations which have smaller ice throughflow and surface slopes. In simulation s1, which has the least amount of ice, maximum sliding speed is less than 100 m a−1 at the Sargans diffluence.

In the lobes, sliding speed generally diminishes quickly as ice radially spreads into the lowlands but sliding speed patterns there also follow topographic patterns of troughs and bedrock highs. For example, sliding speed is relatively high (150 m a−1) in the two main branches (Limmat and Glatt valleys) of the Linth lobe. In the larger Rhine lobe, sliding is significant along the Lake Constance trough and the Thur valley near Frauenfeld (see Fig. 3 for locations), with speeds in excess of 100 m a−1. These patterns are more pronounced for simulation s5, which has the highest ice flux.

Basal shear stress, another indicator of basal condition, is shown in Fig. 10 and is calculated using Eq. (6). For the simulation based on the ice surface reconstruction of (Fig. 10a), basal shear stress is extremely low in the Rhine and Linth lobes, not exceeding 0.02 MPa. Outside of the lobes, values in excess of 0.3 MPa are found in areas with high slope gradients and high surface velocity (see Fig. 7a). In other numerical simulations (s1–s5), basal shear stress depends on the mass balance. For simulation s1 with the smallest mass balance gradient, computed basal shear stress is low in the lobe (up to 0.05 MPa) and up to 0.1 MPa in Alpine valleys. Higher values are obtained for other simulations with higher mass balance gradients and higher ice mass turnover. For these simulations, basal shear stress in the lobe can reach up to 0.1 MPa, twice the value computed for s1, and up to 0.3 MPa in Alpine valleys. The large difference in basal shear stress between lobes and Alpine valleys is consistent among all present simulations and is also a feature of earlier models . In high-elevation areas where ice is frozen to the bed, basal shear stress can also exceed 0.25 MPa. High values of basal shear stress (up to 0.3 MPa) are in accord with results from other numerical modeling studies and with local measurements beneath glaciers . Also, the basal shear stress at the margin of the lobes is high owing to a transition from temperate (sliding) to frozen (no slip) basal conditions.

Figure 10Basal shear stress in megapascals for (a)  (steady-state Stokes solution) and (b–f) simulations s1–s5.

Finally, the ratio of the sliding speed to the surface speed is shown in Fig. 11 as an indicator of internal deformation. This ratio is 0 for the geomorphic reconstruction (Fig. 11a) since no sliding was assumed there. A ratio near 1 indicates little or no internal deformation. The ratio of sliding to surface speed shows a clear pattern for all simulations. Sliding is proportionally greater (0.8 of surface speed) in the lobes than in the main trunks of valleys (0.5–0.7). In the lobe, ice flows almost like a plug with fewer vertical gradients in velocity owing to sliding and to low surface slopes resulting in low driving stresses. The value of the ratio, however, is not homogeneous in the lobe and mimics patterns of basal temperature with less sliding with respect to surface speed at which ice is colder. Sliding ratio increases outward from the center of the lobe up to the cold margin. Unsurprisingly, thicker ice with a temperate, sliding base has a larger component of ice deformation than thinner ice. In the main Alpine valleys, sliding is about 0.5 to 0.6 of surface speed. Slightly higher values are obtained for simulation s1, which has smaller surface gradients and smaller driving forces and thus less internal ice deformation.

Finally, it should be noted that values of sliding speeds, basal shear stresses, and ratios of surface speed to sliding speed are obtained with a linear sliding law and are subject to a cautious interpretation.

Figure 11Ratio of sliding to surface speed for (a)  (no sliding assumed) and (b–f) simulations s1–s5.

4 Discussion

Numerical solutions simulated at least 1500 years of glacier evolution and more than 3000 years in one case (see Table 2). Transient Stokes flow simulations are computer intensive and time consuming, practically limiting the number of runs to be performed and the simulation time. Only two simulations achieved near-steady-state conditions (simulations s1 and s3). For the other three, positive net mass balance (see Table 2) indicates that the ice mass is still increasing and that, at equilibrium, the ice margin would have extended further north than the LGM moraines. This is already the case for simulation s4 (see Fig. 4e). The configurations of the Rhine glacier system for simulations s2, s4, and s5 can be interpreted as the maximum extent at the Rhine glacier at the LGM just prior to its retreat as a result of a global warming of the climate and the termination of the last glaciation.

Although our objective was initially to obtain steady-state configurations of the Rhine glacier under different climate scenarios, the actual Rhine glacier may never have reached equilibrium conditions at the LGM. The climate record prior to the LGM indicates large oscillations as shown by the oxygen isotopic record (δ18O) in Greenland and also from several different proxies in Europe and in the Alps . Two recent speleothems from a cave near Bern, Switzerland, also display these oscillations . These century- to millennia-long oscillations, called Dansgaard–Oeschger (DO) events, occurred repeatedly during Marine Isotope Stage 3 (MIS3, 60–30 ka BP) . Probably as a result of the climate variability during the LGM period, the Rhine glacier left two closely spaced terminal moraines in the Rhine lobe . This twofold maximum is also observed elsewhere in the Alps, on the Tagliamento glacier in the southern Alps in Italy and in Austria . Formation of these two moraines may have been separated in time by less than a couple of thousand years for the Rhine glacier (see Keller and Krayss2005b). Thus, glacier dynamics was likely never in equilibrium with the mass balance. In addition, thermal equilibrium takes significantly more time to establish than dynamic equilibrium. Thus, is it more than likely that the Rhine glacier at the LGM was not a steady-state configuration. This partially validates the use of nonsteady solutions as illustrations of possible Rhine glacier configurations at the LGM.

Our five simulations, although sampling a wide range of ELA and mass balance gradients, do not exhaust the nearly infinite possible combinations of climatic parameters (and glaciological parameters) that could yield other different Rhine glacier configurations. In that sense, the five simulations are part of a continuum of solutions for the Rhine glacier at the LGM. By using a wide range of parameter values for the mass balance gradients and the ELA, we have attempted to bracket as disparate as possible ice configurations of the Rhine glacier at the LGM. Other combinations of climatic and glaciological parameters, however, may result in other ice dynamics configurations not obtainable with the parameters we used. Due to the nature of the model (transient Stokes flow at high resolution for several thousand years), it is not feasible to perform a parametric study and explore more than just a handful of numerical simulations.

Figure 12Comparison of geomorphic reconstruction with extent of temperate ice for simulation s5. (a) Geomorphic ice surface elevation of . Brown represents ice-free peaks above the ice surface; (b) ice surface elevation of simulation s5 with the contour line indicating the position of the temperate–cold transition along mountain flanks and colored according to elevation (same scale as ice surface). (c) Same as (a) with temperate–cold transition of (b) superimposed. All elevations are in meters.

## 4.2 Glacial extent and ice thickness: comparison with geomorphic reconstructions

Two simulations, s2 and s5, have ice cover areas that nearly match the geomorphic reconstruction (16 339 and 15 831 km2, respectively, against 16 400 km2), but are out of equilibrium with ice extent still increasing. For these two simulations, however, the Hörnli ridge is ice covered and the western part of the Rhine lobe is less extensive when compared to the geomorphic reconstruction. The former is due to overall thicker ice for these two simulations. The latter may be due to using the present-day depth of Lake Constance as the basal topography for all of our simulations. Sediments below Lake Constance, however, were partially excavated prior to the LGM as shown by several geologic cross sections based on borehole data . A deeper basal topography along the Lake Constance trough would have channeled more ice along a southeast to northwest path, easing ice advance further to the northwest while reducing the flux of ice in the northeastern part of the Rhine lobe.

In addition to glacial extent, glacial volume may help discriminate among numerical results. Except for simulation s1, all other simulations indicate a thicker glacier (by at least 500 m) than the geomorphic reconstructions. and estimate an ice volume of ∼6500km3, similar to simulation s1 (5500 km3) but significantly less than simulations s2–s5 (10 000 to 15 000 km3; see Table 2). Thicker ice for the Rhine glacier was also obtained by and in a global reconstruction of the ice mass over the European Alps at the LGM, regardless of the precipitation and temperature offsets used. Although simulation s1 yields an ice profile similar to the geomorphic reconstruction (see Fig. 6), there is a discrepancy between simulation s1 mass balance parameters and independent climatic reconstructions for the LGM. For simulation s1, mass balance gradients necessary to obtain such a low-profile Rhine glacier ice surface are very small: 0.1 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}\right)\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$ and 0.025 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{100}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$ for the ablation and accumulation zones, respectively. This value of ablation gradient is identical to the value obtained by using glacier shear stress and mass balance consideration on the reconstructed LGM glacier surface topography along a flow line that extended from the Vorderrhein to Lake Constance. With the ELA set at 1200 m for simulation s1, an ablation gradient of 0.1 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}\right)\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$ yields an ablation rate at the terminus (using 550 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ as a reference) of 0.65 m a−1. The ELA together with the ice surface elevation distribution in the accumulation zone (Fig. 6) yields an average accumulation rate of 0.19 m a−1. The ablation rate can be converted back to a temperature with the model of Reeh (1991). Using a PDD factor of 6 mm PDD−1 and a mean amplitude difference between summer and an annual temperature of 20 C yields a LGM summer temperature at the terminus of about 0 C (Reeh1991, Fig. 2). Higher PDD factors, as perhaps are found in Alpine glaciers (Hock2003), would result in a lower temperature. The LGM accumulation rate is about 10 % to 20 % of today. These numbers are representative of an extremely cold and dry climate.

These numbers seem to be too small and inconsistent with LGM climate as reconstructed from pollen and other proxy records (e.g., Wu et al.2007). Both and report summer temperatures at the elevation of the glacier terminus of 3.2 and 7 C, respectively, which translates, using the PDD model of Reeh (1991) with a PDD factor of 6 mm PDD−1, into 2.1 and 3.9 m a−1 of melt at the terminus. Higher PDD factors would increase these values. Using and values of ELA (see Table 2), these numbers yield ablation gradients significantly higher than those of simulation s1: 0.53 and 0.66 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}\right)\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$ for Benz-Meier and Keller and Krayss, respectively, instead of 0.1 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}\right)\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$ for simulation s1. For the accumulation area, using rough estimates of precipitation of and the area of the accumulation zone based on hypsometric curves yields gradients of 0.034 and 0.06 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}\right)\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$ for Benz-Meier and Keller and Krayss, respectively, instead of 0.025 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}\right)\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$ for simulation s1. Gradients calculated on the basis of or suggest warmer and wetter conditions, closer to values used in simulations s2–s5 (see Table 2). These simulations, however, yielded much thicker ice.

The mass balance ablation gradient of simulation s1, the simulation that best matches the Rhine glacier reconstructed from geomorphic interpretation, is smaller than any present-day values measured in Greenland. Furthermore, the gradients of simulation s1 are only slightly larger than measured in the McMurdo Dry Valleys of Antarctica. The compilation of of surface mass balance gradients in the ablation area of the Greenland Ice Sheet on glaciers and ice caps in Greenland's coastal areas indicates a range between 0.2 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{m}{\right)}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$ in the very north and 0.45 m (100 m) a−1 in the central to southwest as well as the very south of Greenland, about 2 to 4.5 times larger than values obtained in simulation s1. Antarctic mass balance ablation gradients of three Antarctic dry valley glaciers indicate values between 0.02 and 0.05 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{m}{\right)}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{a}}^{-\mathrm{1}}$; the latter is only one-half of the value used in simulation s1. Estimated melt rates at the terminus of the Greenland Ice Sheet and of local glaciers at an altitude similar to the margin of the Rhine glacier lobe range between 1 and 5 m a−1, 1.5 to 7 times larger than the value used in simulation s1. The marginal melt rate in s1 is close to values found in the Antarctic Dry Valleys . Extremely low to even inverse mass balance gradients can locally be produced by extreme shadow (glacier margins in deep-cut valleys), heavy debris cover, or strong snow redistribution with deposition of snow at the ice margin. For the Rhine glacier with a large and predominantly clean piedmont lobe , however, such special conditions are hardly plausible. Thus, although simulation s1 yields an ice thickness comparable with geomorphic reconstructions, its climate may be too extreme in view of our current understanding of climate conditions near the margin of the Rhine glacier at the LGM.

Estimating LGM mass balance gradients and ablation rates based on reconstructed temperature and precipitation patterns results in values close to the ones used in simulations s2 through s5. These simulations, however, yielded much thicker ice, up to 700 m more, than the geomorphic maps of , , or . Based on this comparison alone, our simulation results indicate that either mass balance gradients were high and the Rhine glacier was thicker or gradients were small with the implication that the climate was colder and drier than commonly accepted for the LGM. If the former is correct, this would mean that the trimline used to delineate the glacier extent in the accumulation area, and from which ice surface elevation contour lines are usually drawn , does not represent the ice surface but an englacial cold–temperate ice transition layer that indicates the limit of glacial erosion . In several, more polar paleo-ice sheets (British Isle, Scandinavia, Canada), several studies have pointed out the existence of uneroded bedrock surfaces beneath ice sheets that have survived multiple glaciations under frozen, cold-based basal conditions . The same situation is also believed to exist in some places under Antarctica . From a geomorphic point of view, whether cold ice was present above the trimline during the LGM remains a subject of debate but has not been entirely ruled out (e.g., Wirsig et al.2016). Thus, if instead of the numerically computed ice surface, the cold–temperate ice transition that intersects the basal topography is used for comparing the numerical simulations and the geomorphic reconstruction, the fit between the geomorphic ice surface and the englacial transition is significantly different. Figure 12 shows three images that attempt to compare the geomorphic reconstruction with this cold–temperate ice transition. In Fig. 12a, the ice extent of the geomorphic reconstruction of is shown with visible ice-free peaks above the ice surface (in brown). The boundary of these peaks represents the observed cold–temperate transition (the trimline). Figure 12b shows simulation s5 ice cover (here the ice is so thick that almost no peak protrudes above the ice surface) with the outline of the computed cold–temperate transition colored according to its elevation using the same color scale as the ice surface. Finally, Fig. 12c shows this same cold–temperate transition superimposed on top of the geomorphic reconstruction (Fig. 12a). The elevation of the cold–temperate transition is significantly lower than the ice surface of simulation s5. In the Alps near the foreland, the transition overlaps with the ice-free peaks of the geomorphic map of , which emerge from the ice surface (e.g., Alpstein, Churfirsten; see Fig. 12c). The fit between the trimline of the geomorphic reconstruction (Fig. 12a) and the temperate–cold transition of simulation s5 (Fig. 12b) is nearly perfect. Higher up in the Alps like, for example, along the Glarus Alps, the cold–temperate transition of simulation s5 is at an elevation significantly lower than the trimline of the geomorphic reconstruction (see Fig. 12c) probably because of the accumulation of cold ice at high elevation in simulation s5 that moves englacially. There, the cold–temperate transition underestimates the elevation of the trimline. This discrepancy could easily be explained by the cold temperature at the ice surface used to drive the numerical simulations. Under a warmer and wetter climate that preceded the LGM, thicker temperate ice may have existed in this region with a cold–temperate ice transition located at a higher elevation along the mountain flanks, closer to the measured elevations of the trimlines. Exploring this scenario would require modeling the evolution of the Rhine glacier over longer time periods with a fluctuating climate, something difficult to do with our Stokes flow model. The discrepancy among climate, Rhine glacier ice thickness in the accumulation zone, and trimline reconstruction deserves further attention.

Table 4Annual ice fluxes for the most important glaciers feeding the Linth and Rhine lobes. Cross sections (gates) at which fluxes are computed are shown in Fig. 13a. For the simulation, fluxes are calculated after 50 years of transient simulation with Stokes flow with an evolving ice surface.

## 4.3 Rhine glacier dynamics

### 4.3.1 Surface speed

Our numerical simulations of the Rhine glacier system at the LGM indicate that the ice mass had a relatively slow mass turnover with maximum surface speeds ranging between 180 and 535 m a−1, values that reflect overall small mass balance gradients. Highest flow speeds are obtained for simulation s5, which has the highest rate of accumulation. These surface speeds are significantly smaller than the maximum surface speed estimated from the geomorphic reconstruction. Differences between the geomorphic (Fig. 7a) and numerical reconstructions (Fig. 7b–f) are due to the evolution of the ice-free surface for more than 1000 years in our numerical reconstructions driven by the mass balance forcing (accumulation or ablation) and by ice flow: with time, surface slope gradients become smoother and are glaciologically consistent with the surface mass balance boundary condition, giving rise to the expected patterns of complex ice flow along the main glacial valleys and their tributaries, and with the highest surface speeds in the main valleys in which ice is deepest and at confluences. Maximum surface velocity in simulations s1 to s2 is 3 to 10 times smaller than the results of either or our own reconstruction based on his ice surface elevation. The opposite is true for surface velocity in the lobes. Average computed surface velocities for simulations s1–s5 range from 28 to 79 m a−1, about 30 to 80 times greater than computed for the geomorphic reconstruction using Stokes flow (<1m a−1; see ${\stackrel{\mathrm{‾}}{V}}_{\text{lobe}}$, Table 3). Such low surface speeds for the geomorphic reconstruction are due to (i) no sliding and (ii) a nearly flat lobe. An earlier calculation by that included sliding yielded a value of 25 m a−1.

### 4.3.2 Flow patterns and fluxes

Figure 7 shows surface velocity but also included ice flow lines that delineate major ice-forming basins of the Rhine glacier, from the Ill glacier in Austria in the east, to the Linth glacier that fed the western part of the Linth lobe. To obtain streamlines separating basins, streamlines were either computed forward from a confluence of two glaciers (Linth–Walensee, Rhine–Ill) or backward from a diffluence (Rhine–Walensee at Sargans). The streamlines were generated along a vertical line to follow ice motion at different ice depths. Because ice paths vary with depth, streamlines originating at different ice depths sometimes diverge as, for example, the streamlines backtracked from the Sargans diffluence (in black in Fig. 7).

Results indicate that the Rhine lobe is made up of ice originating from the Valser Rhine, Hinterrhein, Landquart, and Ill glaciers. Ice from the Vorderrhein moves through the Walensee glacier towards the Linth lobe. This differentiation varies somewhat with the thickness of the ice in the accumulation zone: some ice from the Vorderrhein makes it into the Rhine lobe for the thin-ice simulation (s1) while some ice from the Valser Rhine makes it to the Linth lobe for the thickest-ice simulation (s5).

In the Linth lobe, the Linth–Walensee ice divide sometimes passes west of Pfannenstiel (see Figs. 1 and 3 for locations) through the Limmat valley (simulations s1 and s3) and sometimes east of Pfannenstiel through the Glatt valley (simulations s2, s4, and s5). Given the variability, it it possible that this glacier divide was more prone to ice flow switches depending on climatic and glaciological conditions (e.g., local accumulation rates in the Linth and Rhine basins, basal conditions). Also of interest, the Schussen lobe, a minor lobe within the greater Rhine lobe is made up of ice from the Rhine basin and not from the Ill glacier coming from Austria, except for in simulation s4 (Fig. 7e), despite its location in the northeastern part of the Rhine lobe.

Additional insight about flow patterns can be gained by randomly generating tracers in a basin and following their paths downstream. Figure 13 shows ice flow lines colored according to the source region of the ice. Flow lines are calculated from the three-dimensional velocity field. For the geomorphic reconstruction (Fig 13a), the flow lines are based on the velocity field computed after 50 years of transient evolution of the ice surface because the steady-state velocity field (see Fig. 7a) does not yield coherent streamlines.

Figure 13Three-dimensional flow lines for simulations (a)  after 50 years of free surface evolution and (b–f) s1–s5. Colors indicate origin of ice: (green) Linth glacier; (red) Vorderrhein; (yellow) Hinterrhein; (white) Landquart valley; (orange) Ill glacier. Black lines in (a) indicate the four gates through which ice fluxes are calculated (see Table 4). Background color is elevation of basal surface.

The flow lines computed for (geomorphic reconstruction) indicate that ice coming from the eastern glaciers (Ill, Landquart) dominates in the Rhine lobe (Fig. 13a). In contrast, all numerical simulations (s1–s5) indicate that the Rhine lobe consists mostly of ice from the Hinterrhein and the Vorderrhein (Fig. 13b–f). For the geomorphic reconstruction, ice originating from the Hinterrhein has an almost east-to-west flow direction in the Rhine lobe. For the numerical simulations (s1–s5) flow of ice from the Hinterrhein has a southeast-to-northwest orientation. Ice from the Landquart valley is found in the area of Untersee in the geomorphic reconstruction while is it mainly located north towards the Schussen valley in Germany in numerical simulations. Ice from the Ill valley is found as far west as the Überlingen branch of Lake Constance in the geomorphic reconstruction while it is mainly found further east of Schussen valley in the numerical simulations. In the Linth lobe, ice to the west of the Hörnli ridge originates form the Linth glacier while ice to the east in the Glatt valley originates from the Walensee glacier. In simulations s1 and s3 (lowest mass balance gradients of all simulations), some ice from the Walensee flows to the west of the ridge, indicating less ice from the Linth glacier in the eastern part of the model.

Ice thickness also influences flow paths, with thicker ice allowing movement of ice across high passes and ridges in the Alps. While flow is constrained by main valleys in simulation s1 with the lowest mass balance gradients and smallest ice thicknesses, the ice surface in other simulations is significantly higher, allowing ice to move more directly from the accumulation to the ablation across mountain ranges through high passes, for example across the Alpstein and Churfirsten ranges. Despite these differences, however, the overall patterns of ice flow remain strongly constrained by the general complex topography of the Alpine valleys.

Three-dimensional velocities across the different basins can help estimate annual fluxes of ice to the Linth and Rhine lobes (Table 4). Fluxes depend on the rate of accumulation. The range of estimated ice flux can vary by a factor of 5 for the Rhine glacier between the simulations with the smallest accumulation rate (s1, 0.93 km3) and the highest rate (s5, 4.8 km3). The Rhine glacier, because of its large basin area, brings most ice to the Rhine lobe (about 80 %), with little variation among the simulations (s1–s5), and with the remaining ice in the Rhine lobe coming from the Ill glacier. For the Linth lobe, the Walensee glacier brings slightly more ice (0.27 to 1.4 km3, 56–67 % depending on the simulation) than the Linth glacier (0.11 to 1.1 km3). Fluxes computed for the geomorphic reconstruction of yield values close to the driest simulations (s1) for the Linth lobe and closer to intermediate simulations (s2 and s3) for the Rhine lobe. The origin of the ice in the lobe, however, is different, with equal fluxes of ice coming from the Rhine and the Ill glaciers in the Rhine lobe and the Linth and the Walensee glaciers in the Linth lobe. Given the vastly different basin sizes of these glaciers, these numbers do not appear to be realistic. An estimate of ice flux entering the Rhine lobe by yielded a value of 1.86 km3, a value higher than the driest simulation (s1, 1.1 km3) but substantially smaller than the next two wetter simulations (s2 and s3, 3.3 and 3.1 km3, respectively).

### 4.3.3 Lobe dynamics

As ice from the Rhine glacier enters the Alpine foreland through the Lake Constance basin, it fans out in a lobe over 60 km long and 100 km wide. In the lobe, ice flows over undulating topography with numerous troughs and valleys separated by topographic highs (see basal topography in Fig. 3). Sliding speed (Fig. 9) is highest in these troughs in which ice is thicker (Fig. 5) and basal ice is at the melting point (Fig. 8). In the Rhine and Linth lobes, valley troughs with zones of overdeepenings focus ice flow while bedrock highs are zones of lower ice speeds.

Because of topography, sliding speed does not decrease monotonically as ice fans out to the margin in the Rhine lobe. Zones of relatively fast sliding occur near the margin owing to topographic effects and also to ice convergence. In simulations with the thickest ice and with the most extensive glacier cover (s2, s4, and s5), ice velocities, particularly sliding speeds, are sometimes higher near the margin than further up-glacier. This is the case for the Thur valley near the terminus and for Überlingersee (see Figs. 9c, e, f and 1 and 3 for locations). For the Thur valley near the terminus, convergence of ice from Lake Constance and from the upper part of the Thur valley into a low-lying area causes ice to speed up. For the Überlingersee, a narrow trough forces ice to accelerate.

In the Linth lobe the situation is simpler, with ice from the Linth and Walensee glaciers flowing between well constrained ridges that separate the Limmat and Glatt valleys. The Albis and Uetliberg mountains delineate the lobe's western flank. In the center, the Pfannenstiel mountain forces ice to separate and flow around it, and the Hörnli ridge separates the Linth lobe from the Rhine lobe to the east. Except for the value of fluxes and the streamline that separates ice originating from the Linth or Walensee glaciers, flow patterns in the Linth lobe do not change significantly among the different simulations.

Figure 14Extent and thickness in meters (blue to red color map) of temperate basal ice for (a)  after 50 years of free surface evolution and (b–f) simulations s1–s5. Background gray color scheme shows land surface elevation.

## 4.4 Basal conditions

In accord with earlier work , our numerical simulations indicate that the Rhine glacier system is a polythermal glacier. Basal ice is temperate across much of the lobes and in the lower reaches of Alpine valleys and cold elsewhere (see Fig. 8). The temperate basal ice layer, however, is thin, about 100 m in the Alpine valleys and less than 30 m in the lobe (Fig. 14). This represents only a small fraction of the total ice thickness, about 6 % to 7 % in both valleys and the lobes. Volumetrically, cold ice represents roughly 90 % of the total bulk ice but the temperate ice at the bed controls sliding and thus the overall dynamics of the glacier. Contrary to expectations, the temperate basal ice does not grow in thickness down-glacier. The temperate basal ice layer is thickest in the Rhine valley below the Sargans diffluence, probably because of flow constriction associated with the valley width and due to higher viscous strain heating there. Increased flow velocity and strain rates cause more heat to be generated englacially, bringing cold ice near the bed to the melting temperature as noted in earlier modeling studies and in conceptual models . This strain heating effect is shown in Fig. 15, which displays the vertically integrated strain heating projected onto the bed. Significant differences in magnitude exist between the thin-ice simulations (simulation with geomorphic ice surface reconstruction and simulation s1; Fig. 15a, b) and all other simulations with thicker ice (simulations s2–s5; Fig. 15c–f) but patterns of high strain heating are similar.

Figure 15Vertically integrated strain heating in watts per square meter for (a)  after 50 years of free surface evolution and and (b–f) simulations s1–s5.

The Rhine glacier, from the confluence of the Vorderrhein with the Hinterrhein to Lake Constance, is where strain heating is highest, with values in excess of 1 W m−2 for simulations s2–s5. In comparison, geothermal heat flux is of the order of 0.1 W m−2 and does not play a significant role in determining basal temperature. Geothermal heat flux, which varies between 0.06 and 0.12 W m−2 is not correlated with basal temperature: temperate basal conditions are found up-valley in the Alps where the geothermal heat flux is low in comparison to the Swiss lowlands. Flow of ice through Alpine valleys, past constrictions, and around bedrock bumps increases shear deformation and viscous strain heating. Strain heating is larger at the flanks of Alpine valleys than in the center, reflecting higher shear strain rates at which changes in sliding speeds are greater (see Fig. 9). Also, viscous strain heating is large at confluences (e.g., Linth and Walensee, Ill and Rhine) and where ice flows across ridges such as the Appenzell foothills (near where the Rhine glacier enters the Lake Constance basins) or past bedrock bumps like Fläscher Berg near the Sargans diffluence (see Fig. 3 for locations). In contrast, where ice is cold (high Alpine areas) or moves nearly like a plug (Rhine and Linth lobes), viscous strain heating is negligible.

Zones of high viscous strain heating are also zones of higher basal shear stress (see Fig. 10). This is not surprising since areas with high viscous strain heating are the same areas that undergo high shear strain rate at the glacier base. Valley flanks where basal ice temperature transitions from temperate to cold and where there are bedrock ridges and bumps and confluences are zones in which basal shear stress is highest. Bedrock protuberances act like sticky spots (Alley1993), increasing basal resistance. Across these obstacles, basal shear stress can reach values in excess of 0.25 MPa (Fig. 10). Where ice is cold and thus sticks to the base, basal shear stresses are also high. In contrast, where it is mostly temperate like in the Rhine and Linth lobes, basal shear stresses are significantly smaller, with maximum values between 0.05 and 0.1 MPa for thin- and thick-ice simulations, respectively.

Figure 16Effect of basal friction on basal temperature, sliding speed, and basal ice thickness for simulation s5 at t=1540 a. (a–c) No basal friction; (d–f) with basal friction. (a, d) Basal temperature. (b, e) Sliding speed. (c, f) Thickness of temperate basal ice layer. Basal friction was turned on for 200 years. Note also the slightly larger extent of ice for the simulation with basal friction on.

In one simulation we included heat generated by basal friction (see Eq. 8). Figure 16 shows the effects of this additional basal heat source term on basal conditions. Frictional basal heating generates about 0.02 W m−2 on average at the bed of the glacier, increasing basal temperature to the melting temperature everywhere in the lobe, reducing basal resistance due to cold sticky spots. This can explain the increase in the sliding speed by a few tens of meters per year everywhere on the Swiss Plateau, which causes the margin to advance several kilometers further north over the 200 years of the simulation. The effect on the thickness of the temperate ice layer, however, is minimal.

5 Conclusions

Using a fully coupled thermodynamics ice flow model (Elmer/Ice), we investigated the dynamics of the Rhine glacier at the LGM for several climate scenarios characterized by different combinations of ablation and accumulation area mass balance gradients and compared our simulations of ice thickness and flow (where possible) to geomorphic reconstructions and earlier glaciological studies. Our numerical simulations confirmed results of earlier studies but provided more details of the complex three-dimensional ice flow patterns. In all climate scenarios tested, the Rhine glacier is polythermal with temperate basal ice in Alpine valleys downstream of the Vorderrhein–Hinterrhein junction, and in most of the Rhine and Linth lobes. In all simulations, cold and dry climatic conditions yielded a LGM Rhine glacier that moved relatively slowly (maximum ice surface speeds of 150 to 530 m a−1 with even slower speeds in the lobe (average surface speed of 25 to 80 m a−1). Using a linear sliding rule, sliding patterns mimicked surface speeds, with the highest values at the bottom of Alpine valleys. Sliding speed decreased as piedmont glaciers fanned out in the lobes, but topographic lows that follow present-day river drainage focused ice flux and yielded complex patterns of ice flow in the lobes. Basal shear stress in the lobe was smaller than 0.1 MPa for the simulations with high mass balance gradients and down to 0.025 MPa when prescribing gradients with a low mass balance. The latter value is similar to earlier calculations. Higher values, up to 0.3 MPa, were obtained higher up in the glacier and at the margin. In the lobe, a cold margin persisted owing to cold climatic conditions, and ice above several topographic highs was also frozen to the bed. These highs correspond to old uneroded gravel deposits, several hundred to 1 million years old, indicating that parts of the bed resisted glacial erosion. Relatively high flow speeds in some areas near the glacier margin are caused by convergence of ice from different parts of the lobe over troughs such as in the area of the Thur valley.

Our results indicate that the geomorphic ice surface in the accumulation area reconstructed from trimlines can only be matched assuming an extremely cold and dry climate, substantially colder and dryer than estimated for the LGM from proxy records. New pollen records that span the LGM from lakes north of the Swiss Alps (e.g., Duprat-Oualid et al.2017) could help to keep or reject this extreme cold and dry climate scenario needed to numerically reproduce the elevation of the geomorphically reconstructed ice surface at the LGM. Alternatively, simulations that used a less extreme cold–dry climate yielded an ice thickness 500 to 700 m more than the geomorphic reconstruction. This discrepancy can be explained if the trimline and the derived geomorphic ice surface are interpreted as a cold–temperate ice surface above which glacial erosion did not occur, leaving no erosional features on mountain valley flanks. Dating of rock surfaces above the trimline may help resolve this issue.

Our numerical results indicate that transection glaciers such as the Rhine glacier have complex three-dimensional flow patterns that are difficult to reduce to two dimensions. The coupling of ice flow and thermodynamics yields patterns of basal conditions that are inherently three-dimensional, with patches of cold ice within large areas of temperate basal conditions in the lobes. Early studies of the Rhine glacier captured some of its main features but lacked the precision needed to fully understand the coupling among ice flow, climate, and basal thermal conditions. Future models that include the effects of permafrost on basal sliding, and how ice loading and distribution influence patterns of ground water flow, are needed to better understand environmental aspects at the LGM and the characteristics of the subsurface during ice ages in view of the long-term safety of radioactive waste repositories in northern Switzerland under ice house conditions.

Data availability
Data availability.

No new data were acquired to run the simulations.

Author contributions
Author contributions.

DC, WH, and UHF designed the study. DC performed the numerical simulations with the help of FGC. DC carried out the analysis and interpretation of the results with the help of all coauthors. DC wrote the manuscript with input from all coauthors.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

The project was funded by the Swiss Cooperative for the Disposal of Radioactive Waste (Nagra), Wettingen, Switzerland. Most of the computations presented in this paper were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr, last access: 31 July 2018), which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA) and the Equip@meso project (reference ANR-10-EQPX-29-01) of the program Investissements d'Avenir supervised by the French Agence Nationale pour la Recherche. We thank Kevin P. Norton for providing us with a digital map of his results of isostatic uplift at the LGM. Thanks are due to Olivier Gagliardini, Thomas Zwinger, and Martina Schäfer for help with Elmer/Ice and to Chip Hankley and Michael Ruff for their help with the figures. We also thank Shawn Marshall and the anonymous reviewer for their comments that helped improve this paper.

Edited by: G. Hilmar Gudmundsson
Reviewed by: Shawn Marshall and one anonymous referee

References

Adalgeirsdóttir, G., Gudmundsson, G. H., and Björnsson, H.: A regression model for the mass-balance distribution of the Vatnajökull ice cap, Iceland, Ann. Glaciol., 37, 189–193, https://doi.org/10.3189/172756403781815447, 2003. a

Alley, R. B.: In search of ice-stream sticky spots, J. Glaciol., 39, 447–454, https://doi.org/10.3189/S0022143000016336, 1993. a

Anselmetti, F. S., Drescher-Schneider, R., Furrer, H., Graf, H. R., Lowick, S. E., Preusser, F., and Riedi, M. A.: A 180,000 years sedimentation history of a perialpine overdeepened glacial trough (Wehntal, N-Switzerland), Swiss J. Geosci., 103, 345–361, 2010. a

Ballantyne, C. K. and Stone, J. O.: Trimlines, blockfields and the vertical extent of the last ice sheet in southern Ireland, Boreas, 44, 277–287, https://doi.org/10.1111/bor.12109, 2015. 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ürmian) glacial geomorphology: a case study from Oberschwaben (southern Germany), E&G Quaternary Sci. J., 63, 107–129, https://doi.org/10.3285/eg.63.2.01, 2014. a, b, c, d, e

Becker, P., Seguinot, J., Jouvet, G., and Funk, M.: Last Glacial Maximum precipitation pattern in the Alps inferred from glacier modelling, Geogr. Helv., 71, 173–187, https://doi.org/10.5194/gh-71-173-2016, 2016. a, b

Beghin, P., Charbit, S., Dumas, C., Kageyama, M., and Ritz, C.: How might the North American ice sheet influence the northwestern Eurasian climate?, Clim. Past, 11, 1467–1490, https://doi.org/10.5194/cp-11-1467-2015, 2015. a

Benz-Meier, C.: Der würmeiszeitliche Rheingletscher – Maximalstand, Digitale Rekonstruktion, Modellierung und Analyse mit einem Geographischen Informationssystem, Ph.D. thesis, Universität Zürich, 2003. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac, ad, ae, af, ag, ah, ai, aj, ak, al, am, an, ao, ap, aq, ar, as, at, au, av, aw, ax, ay, az, ba, bb, bc, bd, be, bf

Bini, A., Buoncristiani, J.-F., Couterrand, S., Ellwanger, D., Felber, M., Florineth, D., Graf, H. R., Keller, O., Kelly, M., Schlüchter, C., and Schoeneich, P.: Die Schweiz während des letzteiszeitlichen Maximums (LGM), Bundesamt für Landestopografie swisstopo, 2009. a, b, c, d, e

Blatter, H. and Haeberli, W.: Modelling temperature distribution in Alpine glaciers, Ann. Glaciol., 5, 18–22, 1984. a, b, c

Braakhekke, J., Kober, F., and Ivy-Ochs, S.: A GIS-based compilation of the erratic boulders of the Last Glacial Maximum in the Rhine Glacier system, Tech. rep., Nagra Arbeitsbericht NAB 16–60, 2016. a

Brædstrup, C. F., Egholm, D. L., Ugelvig, S. V., and Pedersen, V. K.: Basal shear stress under alpine glaciers: insights from experiments using the iSOSIA and Elmer/Ice models, Earth Surf. Dynam., 4, 159–174, https://doi.org/10.5194/esurf-4-159-2016, 2016. a

Braithwaite, R. J.: Positive degree-day factors for ablation on the Greenland ice sheet studied by energy-balance modelling, J. Glaciol., 41, 153–160, 1995. a, b, c

Briner, J., Miller, G., Davis, P., Bierman, P., and Caffee, M.: Last Glacial Maximum ice sheet dynamics in Arctic Canada inferred from young erratics perched on ancient tors, Quaternary Sci. Rev., 22, 437–444, https://doi.org/10.1016/S0277-3791(03)00003-9, 2003. a

Cohen, D., Hooke, R. L., Iverson, N. R., and Kohler, J.: Sliding of ice past an obstacle at Engabreen, Norway, J. Glaciol., 46, 599–610, 2000. a

Cohen, D., Iverson, N. R., Hooyer, T. S., Fischer, U. H., Jackson, M., and Moore, P. L.: Debris-bed friction of hard-bedded glaciers, J. Geophys. Res.-Earth, 110, F02007, https://doi.org/10.1029/2004JF000228, 2005. a, b

Creyts, T. T., Ferraccioli, F., Bell, R. E., Wolovick, M., Corr, H., Rose, K. C., Frearson, N., Damaske, D., Jordan, T., Braaten, D., and Finn, C.: Freezing of ridges and water networks preserves the Gamburtsev Subglacial Mountains for millions of years, Geophys. Res. Lett., 41, 8114–8122, https://doi.org/10.1002/2014GL061491, 2014. a

Cuffey, K. and Paterson, W.: The physics of glaciers, Academic Press, Amsterdam, 2010. a

Dansgaard, W., Johnsen, S. J., Clausen, H. B., Dahl-Jensen, D., Gundestrup, N. S., Hammer, C. U., Hvidberg, C. S., Steffensen, J. P., Sveinbjornsdottir, A. E., Jouzel, J., and Bond, G.: Evidence for general instability of past climate from a 250-kyr ice-core record, Nature, 364, 218–220, https://doi.org/10.1038/364218a0, 1993. 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, https://doi.org/10.1016/j.quascirev.2011.12.015, 2012. a

Delisle, G., Caspers, G., and Freund, H.: Permafrost in north-central Europe during the Weichselian: how deep?, in: Proceedings of 8th International Conference on Permafrost, Zurich, Switzerland, Belkama Publishers, Lisse, pp. 187–191, 2003. a

Duprat-Oualid, F., Rius, D., Bégeot, C., Magny, M., Millet, L., Wulf, S., and Appelt, O.: Vegetation response to abrupt climate changes in Western Europe from 45 to 14.7k cal a BP: the Bergsee lacustrine record (Black Forest, Germany), J. Quaternary Sci., 32, 1008–1021, https://doi.org/10.1002/jqs.2972, 2017. 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. Sc. Lett., 384, 57–70, 2013. 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

Engelhardt, H. and Kamb, B.: Basal sliding of Ice Stream B, West Antarctica, J. Glaciol., 44, 223–230, https://doi.org/10.3189/S0022143000002562, 1998. a

Fabel, D., Stroeven, A. P., Harbor, J., Kleman, J., Elmore, D., and Fink, D.: Landscape preservation under Fennoscandian ice sheets determined from in situ produced 10Be and 26Al, Earth Planet. Sc. Lett., 201, 397–406, https://doi.org/10.1016/S0012-821X(02)00714-8, 2002. a

Fischer, U. H., Bebiolka, A., Brandefelt, J., Follin, S., Hirschorn, S., Jensen, M., Keller, S., Kennell, L., Näslund, J.-O., Normani, S., Selroos, J.-O., and Vidstrand, P.: Chapter 11 – Radioactive Waste Under Conditions of Future Ice Ages, in: Snow and Ice-Related Hazards, Risks and Disasters, edited by: Shroder, J. F., Haeberli, W., and Whiteman, C., 345–393, Academic Press, Boston, https://doi.org/10.1016/B978-0-12-394849-6.00011-1, 2015. a

Florineth, D.: Surface geometry of the Last Glacial Maximum (LGM) in the Southeastern Swiss Alps (Graubünden) and its paleoclimatological significance, Eiszeitalter u. Gegenwart, 48, 1998. a, b, c, d, e

Florineth, D. and Schlüchter, C.: Reconstructing the Last Glacial Maximum (LGM) ice surface geometry and flowlines in the Central Swiss Alps, Eclogae Geol. Hel., 91, 391–407, 1998. a

Florineth, D. and Schlüchter, C.: Alpine Evidence for Atmospheric Circulation Patterns in Europe during the LGM, Quartenary Res., 54, 295–308, 2000. a

Fountain, A. G., Nylen, T. H., MacClune, K. L., and Dana, G. L.: Glacier mass balances (1993–2001), Taylor Valley, McMurdo Dry Valleys, Antarctica, J. Glaciol., 52, 451–462, 2006. a, b

Gagliardini, O., Cohen, D., Råback, P., and Zwinger, T.: Finite-element modeling of subglacial cavities and related friction law, J. Geophys. Res., 112, F02027, https://doi.org/10.1029/2006JF000576, 2007. 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, https://doi.org/10.5194/gmd-6-1299-2013, 2013. a, b, c

Genty, D., Blamart, D., Ouahdi, R., Gilmour, M., Baker, A., Jouzel, J., and Van-Exter, S.: Precise dating of Dansgaard-Oeschger climate oscillations in western Europe from stalagmite data, Nature, 421, 833–837, https://doi.org/10.1038/nature01391, 2003. a

Geuzaine, C. and Remacle, J.-F.: Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, Int. J. Numer. Meth. Eng., 79, 1309–1331, 2009. a

Greve, R. and Blatter, H.: Dynamics of ice sheets and glaciers, Springer Science & Business Media, 2009. a

Gross, G., Kerschner, H., and Patzelt, G.: Methodische Untersuchungen über die Schneegrenze in alpinen Gletschergebieten, Zeitschrift für Gletscherkunde und Glazialgeologie, 12, 223–251, 1977. a

Guiot, J., de Beaulieu, J., Cheddadi, R., David, F., Ponel, P., and Reille, M.: The climate in Western Europe during the last Glacial/Interglacial cycle derived from pollen and insect remains, Palaeogeogr. Palaeocl., 103, 73–93, https://doi.org/10.1016/0031-0182(93)90053-L, 1993. a

Haeberli, W.: Permafrost glacier relationships in the Swiss Alps – today and in the past, in: Permafrost Fourth International Conference Proceedings, pp. 415–420, Washington DC, National Academy Press, Fairbanks, Alaska, 1983. a

Haeberli, W.: Zur Glaziologie der letzteiszeitlichen Alpenvergletscherung, in: Klimageschichtliche Probleme der letzten 130'000 Jahre, edited by: Frenzel, B., pp. 409–419, Gustav Fischer Verlag, Stuttgart-New York, 1991. a, b

Haeberli, W.: Glaciological conditions in Northern Switzerland during recent Ice Ages, Tech. rep., Nagra Arbeitsbericht NAB 10-18, 2010. a

Haeberli, W. and Penz, U.: An attempt to reconstruct glaciological and climatological characteristics of 18 ka BP Ice Age glaciers in and around the Swiss Alps, Zeitschrift für Gletscherkunde und Glacialgeologie, 21, 351–361, 1985. a, b, c, d, e, f, g, h, i, j, k, l

Haeberli, W. and Schlüchter, C.: Geological evidence to constrain modelling of the Late Pleistocene Rhonegletscher (Switzerland), in: The Physical basis of Ice Sheet Modelling, 1987. a, b

Haeberli, W., Rellstab, W., and Harrison, W. D.: Geothermal effects of 18 ka BP ice conditions in the Swiss Plateau, Ann. Glaciol., 5, 56–60, 1984. a

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

Hofer, D., Raible, C. C., Merz, N., Dehnert, A., and Kuhlemann, J.: Simulated winter circulation types in the North Atlantic and European region for preindustrial and glacial conditions, Geophys. Res. Lett., 39, L15805, https://doi.org/10.1029/2012GL052296, 2012. a, b, c

Iverson, N. R.: Shear resistance and continuity of subglacial till: hydrology rules, J. Glaciol., 56, 1104–1114, https://doi.org/10.3189/002214311796406220, 2010. a

Iverson, N. R., Cohen, D., Hooyer, T. S., Fischer, U. H., Jackson, M., Moore, P. L., Lappegard, G., and Kohler, J.: Effects of Basal Debris on Glacier Flow, Science, 301, 81–84, https://doi.org/10.1126/science.1083086, 2003. a

Ivy-Ochs, S.: Glacier variations in the European Alps at the end of the last glaciation, Cuadernos de investigación geográfica, 41, 295–315, 2015. a

Jäckli, H.: Die Vergletscherung der Schweiz im Würm maximum, Eclogae Geol. Helv, 55, 285–594, 1962. a, b

Jäckli, H.: Die Schweiz zur letzten Eiszeit, in: Atlas der Schweiz, Eidg. Landestopographie, Wabern-Bern, 1970. a, b

Jóhannesson, T., Sigurdsson, O., Laumann, T., and Kennett, M.: Degree-day glacier mass-balance modelling with applications to glaciers in Iceland, Norway and Greenland, J. Glaciol., 41, 345–358, https://doi.org/10.3189/S0022143000016221, 1995. a

Jouvet, G., Seguinot, J., Ivy-Ochs, S., and Funk, M.: Modelling the diversion of erratic boulders by the Valais Glacier during the last glacial maximum, J. Glaciol., 63, 487–498, 2017. a, b

Keller, O.: Ältere spätwürmzeitliche Gletschervorstöße und Zerfall des Eisstromnetzes in den Nördlichen Rhein-Alpen (Weißbad-Stadium, Bühl-Stadium), vol. 27, Univ. Zürich-Irchel, Geographisches Inst., 1988. a

Keller, O. and Krayss, E.: Die hochwürmzeitlichen Rückzugsphasen des Rhein-Vorlandgletschers und der erste alpine Eisrandkomplex im Spätglazial, Geogr. Helv., 42, 169–178, https://doi.org/10.5194/gh-42-169-1987, 1987. a, b

Keller, O. and Krayss, E.: Datenlage und Modell einer Rhein-Linth-Vorlandvergletscherung zwischen Eem-Interglazial und Hochwürm, GeoArchaeoRhein, 2, 121–138, 1989. a, b

Keller, O. and Krayss, E.: The Rhine-Linth glacier in the upper Würm: A model of the last Alpine glaciation, Quaternary Int., 18, 15–27, 1993. a, b

Keller, O. and Krayss, E.: Die Bodensee-Vorlandvereisung des Rheingletschers im Konstanz-Stadium der letzten Eiszeit, Berichte St-Gallischen naturwissenschaftlichen Ges., 87, 31–40, 1994. a, b

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, b

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, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s

Kelly, J., Buoncristiani, C., and Schlüchter, C.: A reconstruction of the last glacial maximum (LGM) ice surface geometry in the western Swiss Alps and contiguous Alpine regions in Italy and France, Eclogae Geol. Hel., 97, 57–75, 2004. a, b, c

Kirchner, N., Ahlkrona, J., Gowan, E., Lötstedt, P., Lea, J., Noormets, R., von Sydow, L., Dowdeswell, J., and Benham, T.: Shallow ice approximation, second order shallow ice approximation, and full Stokes models: A discussion of their roles in palaeo-ice sheet modelling and development, Quaternary Sci. Rev., 135, 103–114, https://doi.org/10.1016/j.quascirev.2016.01.013, 2016. a, b

Krabbendam, M.: Sliding of temperate basal ice on a rough, hard bed: creep mechanisms, pressure melting, and implications for ice streaming, The Cryosphere, 10, 1915–1932, https://doi.org/10.5194/tc-10-1915-2016, 2016. a

Levermann, A. and Winkelmann, R.: A simple equation for the melt elevation feedback of ice sheets, The Cryosphere, 10, 1799–1807, https://doi.org/10.5194/tc-10-1799-2016, 2016. a

Lindgren, A., Hugelius, G., Kuhry, P., Christensen, T. R., and Vandenberghe, J.: GIS-based maps and area estimates of Northern Hemisphere permafrost extent during the Last Glacial Maximum, Permafrost Periglac., 27, 6–16, https://doi.org/10.1002/ppp.1851, 2016. a

Linsbauer, A., Paul, F., and Haeberli, W.: Modeling glacier thickness distribution and bed topography over entire mountain ranges with GlabTop: Application of a fast and robust approach, J. Geophys. Res.-Earth, 117, F03007, https://doi.org/10.1029/2011JF002313, 2012. a

Löfverström, M., Caballero, R., Nilsson, J., and Kleman, J.: Evolution of the large-scale atmospheric circulation in response to changing ice sheets over the last glacial cycle, Clim. Past, 10, 1453–1471, https://doi.org/10.5194/cp-10-1453-2014, 2014. a

Luetscher, M., Boch, R., Sodemann, H., Spötl, C., Cheng, H., Edwards, R. L., Frisia, S., Hof, F., and Müller, W.: North Atlantic storm track changes during the Last Glacial Maximum recorded by Alpine speleothems, Nat. Commun., 6, 6344, https://doi.org/10.1038/ncomms7344, 2015. a, b, c, d

Lüthi, M. P., Ryser, C., Andrews, L. C., Catania, G. A., Funk, M., Hawley, R. L., Hoffman, M. J., and Neumann, T. A.: Heat sources within the Greenland Ice Sheet: dissipation, temperate paleo-firn and cryo-hydrologic warming, The Cryosphere, 9, 245–253, https://doi.org/10.5194/tc-9-245-2015, 2015. a

Machguth, H. and Cohen, D.: Estimating the surface mass balance gradient of the Rhine Glacier at the Last Glacial Maximum, Tech. rep., Nagra, unpublished. a, b

Machguth, H., Thomsen, H. H., Weidick, A., Abermann, J., Ahlström, A. P., Andersen, M. L., Andersen, S. B., Björk, A. A., Box, J. E., Braithwaite, R. J., Bøggild, C. E., Citterio, M., Clement, P., Colgan, W., Fausto, R. S., Gleie, K., Hasholt, B., Hynek, B., Knudsen, N. T., Larsen, S. H., Mernild, S., Oerlemans, J., Oerter, H., Olesen, O. B., Smeets, C. J. P. P., Steffen, K., Stober, M., Sugiyama, S., van As, D., van den Broeke, M. R., and van de Wal, R. S.: Greenland surface mass-balance observations from the ice-sheet ablation area and local glaciers, J. Glaciol., 62, 861–887, https://doi.org/10.1017/jog.2016.75, 2016. a, b

Marshall, S. J., James, T. S., and Clarke, G. K.: North American Ice Sheet reconstructions at the Last Glacial Maximum, Quaternary Sci. Rev., 21, 175–192, 2002. a

Marshall, S. J., Sharp, M. J., Burgess, D. O., and Anslow, F. S.: Near-surface-temperature lapse rates on the Prince of Wales Icefield, Ellesmere Island, Canada: implications for regional downscaling of temperature, Int. J. Climatol., 27, 385–398, https://doi.org/10.1002/joc.1396, 2007. a

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

Monegato, G., Ravazzi, C., Donegana, M., Pini, R., Calderoni, G., and Wick, L.: Evidence of a two-fold glacial advance during the last glacial maximum in the Tagliamento end moraine system (eastern Alps), Quaternary Res., 68, 284–302, 2007. a

North Greenland Ice Core Project members: High-resolution record of Northern Hemisphere climate extending into the last interglacial period, Nature, 431, 147–151, https://doi.org/10.1038/nature02805, 2004. a

Norton, K. P. and Hampel, A.: Postglacial rebound promotes glacial re-advances – a case study from the European Alps, Terra Nova, 22, 297–302, https://doi.org/10.1111/j.1365-3121.2010.00946.x, 2010. a

Paul, F. and Linsbauer, A.: Modeling of glacier bed topography from glacier outlines, central branch lines, and a DEM, Int. J. Geogr. Inf. Sci., 26, 1173–1190, 2012. a

Penck, A. and Brückner, E.: Die Alpen im Eiszeitalter, Tauchnitz, Leipzig, 1909. a

Pietsch, J. and Jordan, P.: Digitales Höhenmodell Basis Quartär der Nordschweiz – Version 2014 (SGT E2) und ausgewählte Auswertungen, Tech. rep., Nagra Arbeitsbericht NAB 14-02, 2014. a

Pohjola, V. A. and Hedfors, J.: Studying the effects of strain heating on glacial flow within outlet glaciers from the Heimefrontfjella Range, Dronning Maud Land, Antarctica, Ann. Glaciol., 37, 134–142, 2003. 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, 282–305, https://doi.org/10.3285/eg.60.2-3.06, 2011. a, b

Railsback, L. B., Gibbard, P. L., Head, M. J., Voarintsoa, N. R. G., and Toucanne, S.: An optimized scheme of lettered marine isotope substages for the last 1.0 million years, and the climatostratigraphic nature of isotope stages and substages, Quaternary Sci. Rev., 111, 94–106, https://doi.org/10.1016/j.quascirev.2015.01.012, 2015. a

Reeh, N.: Parameterization of melt rate and surface temperature on the Greenland ice sheet, Polarforschung, 59, 113–128, 1991. a, b, c, d

Robin, G. d. Q.: Ice movement and temperature distribution in glaciers and ice sheets, J. Glaciol., 2, 523–532, 1955. a, b

Ryser, C., Lüthi, M. P., Andrews, L. C., Catania, G. A., Funk, M., Hawley, R., Hoffman, M., and Neumann, T. A.: Caterpillar-like ice motion in the ablation zone of the Greenland Ice Sheet, J. Geophys. Res.-Earth, 119, 2258–2271, https://doi.org/10.1002/2013JF003067, 2014. 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 Chronology Part I: Europe, edited by: Ehlers, J. and Gibbard, P., vol. 2, Part 1 of Developments in Quaternary Sciences, 413–418, Elsevier, https://doi.org/10.1016/S1571-0866(04)80092-7, 2004. a

Schoof, C.: The effect of cavitation on glacier sliding, Proc. R. Soc. A, 461, 609–627, https://doi.org/10.1098/rspa.2004.1350, 2005. a

Seddik, H., Greve, R., Zwinger, T., Gillet-Chaulet, F., and Gagliardini, O.: Simulations of the Greenland ice sheet 100 years into the future with the full Stokes model Elmer/Ice, J. Glaciol., 58, 427–440, 2012. a, b, c

Seguinot, J.: Numerical modelling of the Cordilleran ice sheet, Ph.D. thesis, Stockholm University, 2014. a

Seguinot, J., Rogozhina, I., Stroeven, A. P., Margold, M., and Kleman, J.: Numerical simulations of the Cordilleran ice sheet through the last glacial cycle, The Cryosphere, 10, 639–664, https://doi.org/10.5194/tc-10-639-2016, 2016. a

Seguinot, J., Jouvet, G., Huss, M., Funk, M., Ivy-Ochs, S., and Preusser, F.: Modelling last glacial cycle ice dynamics in the Alps, The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-8, in review, 2018. a, b

Shreve, R.: Glacier sliding at subfreezing temperatures, J. Glaciol., 30, 341–347, 1984. a

Speck, C. K.: Änderung des Grundwasserregimes unter dem Einfluss von Gletschern und Permafrost, Ph.D. thesis, Technischen Hochschule Zürich, 1994. a

Spötl, C. and Mangini, A.: U/Th age constraints on the absence of ice in the central Inn Valley (eastern Alps, Austria) during Marine Isotope Stages 5c to 5a, Quaternary Res., 66, 167–175, https://doi.org/10.1016/j.yqres.2006.03.002, 2006.  a

Stokes, C. R., Tarasov, L., Blomdin, R., Cronin, T. M., Fisher, T. G., Gyllencreutz, R., Hättestrand, C., Heyman, J., Hindmarsh, R. C., Hughes, A. L., Jakobsson, M., Kirchner, N., Livingstone, S. J., Margold, M., Murton, J. B., Noormets, R., Peltier, W. R., Peteet, D. M., Piper, D. J., Preusser, F., Renssen, H., Roberts, D. H., Roche, D. M., Saint-Ange, F., Stroeven, A. P., and Teller, J. T.: On the reconstruction of palaeo-ice sheets: Recent advances and future challenges, Quaternary Sci. Rev., 125, 15–49, https://doi.org/10.1016/j.quascirev.2015.07.016, 2015. a

Tarasov, L. and Peltier, W. R.: Impact of thermomechanical ice sheet coupling on a model of the 100 kyr ice age cycle, J. Geophys. Res.-Atmos., 104, 9517–9545, https://doi.org/10.1029/1998JD200120, 1999. a

Tarasov, L., Dyke, A. S., Neal, R. M., and Peltier, W.: A data-calibrated distribution of deglacial chronologies for the North American ice complex from glaciological modeling, Earth Planet. Sc. Lett., 315–316, 30–40, https://doi.org/10.1016/j.epsl.2011.09.010, 2012. a

Thoma, M., Grosfeld, K., Barbi, D., Determann, J., Goeller, S., Mayer, C., and Pattyn, F.: RIMBAY – a multi-approximation 3D ice-dynamics model for comprehensive applications: model description and examples, Geosci. Model Dev., 7, 1–21, https://doi.org/10.5194/gmd-7-1-2014, 2014. a

van den Broeke, M., Bus, C., Ettema, J., and Smeets, P.: Temperature thresholds for degree-day modelling of Greenland ice sheet melt rates, Geophys. Res. Lett., 37, l18501, https://doi.org/10.1029/2010GL044123, 2010. a

van Husen, D.: LGM and Late-glacial fluctuations in the Eastern Alps, Quaternary Int., 38, 109–118, 1997. a

Weertman, J.: Stability of ice-age ice sheets, J. Geophys. Res., 66, 3783–3792, https://doi.org/10.1029/JZ066i011p03783, 1961. a

Wirsig, C., Zasadni, J., Christl, M., Akçar, N., and Ivy-Ochs, S.: Dating the onset of LGM ice surface lowering in the High Alps, Quaternary Sci. Rev., 143, 37–50, https://doi.org/10.1016/j.quascirev.2016.05.001, 2016. a, b

Wu, H. B., Guiot, J. L., Brewer, S., and Guo, Z. T.: Climatic changes in Eurasia and Africa at the last glacial maximum and mid-Holocene: reconstruction from pollen data using inverse vegetation modelling, Clim. Dynam., 29, 211–229, 2007. a

Zoet, L. K. and Iverson, N. R.: Rate-weakening drag during glacier sliding, J. Geophys. Res.-Earth, 121, 1206–1217, https://doi.org/10.1002/2016JF003909, 2016. a, b