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

At the Last Glacial Maximum (LGM), the Rhine glacier in the Swiss Alps covered an area of about 16 000 km². 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}{}^{\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⁻¹ in the lobes and 50 to 250 m a⁻¹ 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.2 MPa). 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⁻². Basal friction acted as a heat source at the bed of about 0.02 W m⁻², 4 to 6 times less than the geothermal heat flow which is locally high (up to 0.12 W m⁻²). 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 1 Map 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 (Preusser et al., 2011) 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 (e.g., Penck and Brückner, 1909; Schlüchter, 1988, 2004; Keller, 1988; Beckenbach et al., 2014). Erratics, tills, and other sediment deposits provided information about flow paths and provenances (e.g., Haeberli and Schlüchter, 1987; Florineth, 1998; Anselmetti et al., 2010; Ellwanger et al., 2011; Braakhekke et al., 2016). 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 (Florineth, 1998; Beckenbach et al., 2014). In the accumulation area, trimlines represent the upper limit of glacial erosion and define the minimum elevation of the ice surface (Florineth and Schlüchter, 1998; Kelly et al., 2004). 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., Jäckli, 1962, 1970; Keller and Krayss, 1987, 1989, 1993, 1994; Benz-Meier, 2003; Kelly et al., 2004; Bini et al., 2009).

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 (Blatter and Haeberli, 1984; Haeberli and Penz, 1985) 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 Penz, 1985). 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 (Florineth, 1998; Florineth and Schlüchter, 2000; Hofer et al., 2012) implied that most of the moisture feeding glaciers in the Alps had a southern, Mediterranean origin (Luetscher et al., 2015). Numerical models of LGM paleo-climate (e.g., Hofer et al., 2012; Löfverström et al., 2014; Beghin et al., 2015) 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⁻¹ for the Rhine lobe (Haeberli and Penz, 1985). Average ice flow velocity through the main outlets of the Rhine glacier was estimated to be less than 200 m a⁻¹ (Keller and Krayss, 2005b), 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 (Haeberli et al., 1984; Delisle et al., 2003; Lindgren et al., 2016). Subsurface temperatures and groundwater flow conditions must have been strongly influenced by the presence of extended surface and subsurface ice (Speck, 1994). 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. 2; Preusser 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 (Blatter and Haeberli, 1984; Haeberli and Penz, 1985) 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-Meier, 2003) 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 (Kirchner et al., 2016), and ignored glacier thermodynamics. An estimate of potential erosion for the Rhine glacier using this approximation (Dürst Stucki and Schlunegger, 2013) 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 (e.g., Benz-Meier, 2003; Dürst Stucki and Schlunegger, 2013) should abide by the fundamental principles of glacier dynamics (Tarasov et al., 2012; Stokes et al., 2015). 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., Haeberli, 1991). 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 (Fischer et al., 2015). 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 (Gagliardini et al., 2013), 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 2 Depth of overdeepened valleys in the Rhine glacier basin (Pietsch and Jordan, 2014). Outlines of maximum glacial extent at the LGM (Bini et al., 2009) 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 (Gagliardini et al., 2013). 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 (Ryser et al., 2014; Kirchner et al., 2016). 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 (Linsbauer et al., 2012; Paul and Linsbauer, 2012) and lake bathymetry where available (Lake Constance and Lake Zurich; Benz-Meier, 2003) and by depressing the topography to reproduce isostatic adjustments at the LGM (up to 130 m) using the model of Norton and Hampel (2010). Figure 3 shows the basal map used for the model.

Figure 3 Basal topography (from Benz-Meier, 2003) 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 Benz-Meier (2003). 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 Benz-Meier (2003) or an ice surface obtained from a previous numerical simulation. The reconstruction of Benz-Meier (2003) is based on the earlier works of Jäckli (1962, 1970) and also of Keller and Krayss (1987, 1989, 1993, 1994). We refer to it as Benz-Meier (2003) despite its earlier origin. A two-dimensional, fixed, unstructured triangular mesh over the domain is created using Gmsh (Geuzaine and Remacle, 2009) 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 \mathit{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}}\mathit{v}+\left[\mathrm{\nabla }\mathit{v}+{\left(\mathrm{\nabla }\mathit{v}\right)}^T\right]\mathrm{\nabla }\mathit{\eta }+\mathit{\rho }\mathit{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}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(EA\left(T\right)\right)}^{-\mathrm{1}/n}{d}^{-(\mathrm{1}-\mathrm{1}/n)},d>d_o,\\ \\ {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(EA\left(T\right)\right)}^{-\mathrm{1}/n}{d}_o^{-(\mathrm{1}-\mathrm{1}/n)},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, d_o is the critical strain rate, $D=\frac{\mathrm{1}}{\mathrm{2}}\left[\left(\mathrm{\nabla }\mathit{v}\right)+{\left(\mathrm{\nabla }\mathit{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)\exp \left({\displaystyle \frac{-Q\left(T\right)}{RT}}\right),\end{array}$$

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

2.3 The heat equation

The heat equation for a viscous fluid is

$$\begin{array}{}\text{(5)}& \mathit{\rho }c\left(T\right)\left({\displaystyle \frac{\partial T}{\partial t}}+\mathit{v}\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 1 Model 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⁻² 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 Medici and Rybach (1995) 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 Medici and Rybach (1995); 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., Shreve, 1984) 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)}& {\mathit{\tau }}_{\mathrm{b}}=C\Vert {\mathit{v}}_{\mathrm{s}}{\Vert }^{m-\mathrm{1}}{\mathit{v}}_{\mathrm{s}},\end{array}$$

where v_s 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; Iverson, 2010), flow over bedrock protuberances with ice–bed separation and water-filled cavities (e.g., Schoof, 2005; Gagliardini et al., 2007; Zoet and Iverson, 2016), or fast ice sheet flows (e.g., Engelhardt and Kamb, 1998). Given the large uncertainties in predicting and modeling glacial sliding (e.g., Marshall et al., 2002; Cohen et al., 2005; Zoet and Iverson, 2016), 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, v_s=0. To simulate this behavior we make the sliding coefficient C in the sliding rule above a function of temperature T (Greve and Blatter, 2009),

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

where C_o and C₁ are the sliding coefficients for temperate and cold conditions, respectively, and $T^{\prime }=T_{\text{pmp}}-T$, where T_pmp is the temperature at the pressure melting point. The parameter γ, sometimes called the sub-melt sliding parameter (Seddik et al., 2012), 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 (Seddik et al., 2012)

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

where n is the normal to the glacier bed, q_geo 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 $\mathit{t}\cdot \mathit{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)}& {\displaystyle \frac{\partial \mathit{\zeta }}{\partial t}}=\dot{b}+w-u{\displaystyle \frac{\partial \mathit{\zeta }}{\partial x}}-v{\displaystyle \frac{\partial \mathit{\zeta }}{\partial y}},\end{array}$$

where ζ is the glacier surface elevation and $\dot{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 ($\dot{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)}& \dot{b}=\left\{\begin{array}{lr}min\left({\dot{b}}_{\text{acc}}^{\text{max}},{\mathit{\beta }}_{\text{acc}}(z-z_{\text{ela}})\right),& z\ge z_{\text{ela}},\\ {\mathit{\beta }}_{\text{abl}}(z-z_{\text{ela}}),& 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 ${\dot{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 (e.g., Becker et al., 2016; Seguinot et al., 2016, 2018; Jouvet et al., 2017). Furthermore, PDD factors needed to compute surface melt rates are known to vary substantially (see Braithwaite, 1995; Hock, 2003; van den Broeke et al., 2010) 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 (Florineth, 1998; Hofer et al., 2012; Luetscher et al., 2015), 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 ${\dot{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 (Haeberli and Penz, 1985) or present-day Greenland values (Machguth et al., 2016) 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 (e.g., Tarasov and Peltier, 1999; Seddik et al., 2012; Seguinot, 2014; Thoma et al., 2014; Jouvet et al., 2017), 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 −12 ^∘C, a value within a range of estimations (−15 to −10 ^∘C; Haeberli, 1983, 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 T_ela and z_ela are the temperature and elevation at the equilibrium line, respectively. In all model simulations ${\mathit{\gamma }}_a=-\mathrm{6}$ ${}^{\circ }\mathrm{C}{\mathrm{km}}^{-\mathrm{1}}$, an intermediate value based on an estimate by Keller and Krayss (2005b) (−7 ${}^{\circ }\mathrm{C}{\mathrm{km}}^{-\mathrm{1}}$) for the Rhine glacier at the LGM and contemporary lapse rates for polar regions (−5 to −4 ${}^{\circ }\mathrm{C}{\mathrm{km}}^{-\mathrm{1}}$; Marshall et al., 2007; Machguth and Cohen, unpublished). Estimates of the equilibrium line altitude at the LGM vary from 944 m (Benz-Meier, 2003) to 1200 m (Haeberli and Penz, 1985). 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 Keller and Krayss (2005b) 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 Robin (1955) (see Cuffey and Paterson, 2010, 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=-\dot{b}z/H$ is, for the accumulation zone ($\dot{b}>\mathrm{0}$),

$$\begin{array}{}\text{(12)}& T=T_{\text{surf}}+z_{\star }{\displaystyle \frac{\sqrt{\mathit{\pi }}}{\mathrm{2}}}{\left[{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\right]}_{\text{bed}}\left[\text{erf}(z/z_{\star })-\text{erf}(H/z_{\star })\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/\dot{b}$. For ablation ($\dot{b}<\mathrm{0}$), the temperature distribution is

$$\begin{array}{ll}{\displaystyle }T=T_{\text{surf}}& {\displaystyle }+{\left[{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\right]}_{\text{bed}}\left\{z_{\star }\exp (z^{\mathrm{2}}/z_{\star }^{\mathrm{2}})\mathcal{D}(z/z_{\star })\right.\\ \text{(13)}& {\displaystyle }& {\displaystyle }\left.-z_{\star }\exp (H^{\mathrm{2}}/z_{\star }^{\mathrm{2}})\mathcal{D}(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 Robin (1955) 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 2 Summary of numerical (this study) and selected geomorphic reconstructions. z_ela is the elevation of the equilibrium line. β_abl and β_acc are the ablation and accumulation mass balance gradients, respectively. ${\dot{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. ${\dot{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. t_s 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 Benz-Meier (2003) for s1, a simulation with β_acc=0.1 $\mathrm{m}(\mathrm{100}\mathrm{m}{)}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$, β_abl=0.2 $\mathrm{m}(\mathrm{100}\mathrm{m}{)}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$, and z_ELA=1200 m that ran for 440 years using the Benz-Meier (2003) reconstruction as the initial condition for s2, and a simulation with β_acc=0.05 $\mathrm{m}(\mathrm{100}\mathrm{m}{)}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$, β_abl=0.2 $\mathrm{m}(\mathrm{100}\mathrm{m}{)}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$, and z_ELA=900 m that ran for 907 years using Benz-Meier (2003) as the initial condition for s3, s4, and s5.

^aBenz-Meier (2003). ^bKeller and Krayss (2005b). ^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 (Haeberli and Penz, 1985; Benz-Meier, 2003; Keller and Krayss, 2005a) 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 (Luetscher et al., 2015) 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 Adalgeirsdóttir et al. (2003) 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: Benz-Meier (2003) and Keller and Krayss (2005b).

Input and output quantities are different for the numerical and the geomorphic reconstructions. For the numerical simulations, model input parameters are the ELA (z_ela), the mass balance gradients (β_abl and β_acc), and the maximum rate of accumulation (${\dot{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 (${\dot{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, t_s, 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, ${\dot{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 ${\dot{b}}_{\text{term}}$ is calculated from values of summer and mean annual LGM temperatures near the 500 m elevation level indicated in both Benz-Meier (2003) and Keller and Krayss (2005b). 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⁻¹ (Braithwaite, 1995), a value used for melting ice in Greenland. Higher PDD values (6–8 mm PDD⁻¹) have been measured on mountain glacier (e.g., Braithwaite, 1995; Hock, 2003; Jóhannesson et al., 1995) but LGM climate conditions in northern Switzerland may have likely been more similar to present-day Arctic conditions. Both Benz-Meier (2003) and Keller and Krayss (2005b) 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 Keller and Krayss (2005b) (0.50 m a⁻¹, Table 2) is higher than the value cited in Keller and Krayss (2005b) (0.30 m a⁻¹) obtained from numbers cited in Haeberli (1991) that assumed LGM precipitation was equal to 20 % of today's value (1.5 m a⁻¹).

2.9 Numerical solutions

The open-source software Elmer/Ice (Gagliardini et al., 2013) 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 4 Ice surface elevation in meters above sea level for (a) Benz-Meier (2003) 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 Benz-Meier (2003) 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 ${\dot{b}}_{\text{net}}=\mathrm{0.02}\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$; see Table 2) and simulation s3 seems close to it (${\dot{b}}_{\text{net}}=\mathrm{0.19}\mathrm{m}{\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 (Weertman, 1961; Levermann and Winkelmann, 2016). 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 km² while simulation s1 is only 12 389 km² (the geomorphic reconstructions yielded values around 16 000 km²). The ice surface area of simulation s5 is similar to the geomorphic reconstruction (16 400 km²). Simulation s4 is slightly larger (17 706 km²); s2 is slightly smaller (15 813 km²). 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 Benz-Meier (2003) (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 Benz-Meier (2003) 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 5 Ice thickness in meters for (a) Benz-Meier (2003) 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 Benz-Meier (2003) and Keller and Krayss (2005b) 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 6 Glacier 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 7 Ice surface speed in meters per year for (a) Benz-Meier (2003) (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⁻¹. 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 Benz-Meier (2003) 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 Benz-Meier (2003). Although our Stokes flow model is different from the simple one-dimensional model of Benz-Meier (2003), our patterns of ice flow are comparable to his reconstruction (see Benz-Meier, 2003, Fig. 6.10) and ice velocities are comparable in magnitude for the most part (e.g., 2930 and 3130 m a⁻¹ 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⁻¹, far exceeding velocities predicted by Benz-Meier (∼2500 m a⁻¹). Table 3 summarizes maximum speeds for all simulations as well as other dynamical quantities described in the following sections.

Benz-Meier (2003)

Table 3 Summary of ice dynamics quantities. V^max 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 Benz-Meier (2003). 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 (Beckenbach et al., 2014) and to thinner ice cover there (see Fig. 5).

Figure 8 Basal temperature relative to the melting point ($\mathrm{d}T=T-T_{\text{pmp}}$) in Kelvin for (a) Benz-Meier (2003) (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 Benz-Meier (2003), 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 9 Ice sliding speed in meters per year for (a) Benz-Meier (2003) (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⁻¹. 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⁻¹ 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⁻¹) 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⁻¹. 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 Benz-Meier (2003) (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 (Haeberli and Penz, 1985). 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 (Brædstrup et al., 2016) and with local measurements beneath glaciers (Cohen et al., 2000, 2005; Iverson et al., 2003). 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 10 Basal shear stress in megapascals for (a) Benz-Meier (2003) (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 11 Ratio of sliding to surface speed for (a) Benz-Meier (2003) (no sliding assumed) and (b–f) simulations s1–s5.

4 Discussion

4.1 Steady state

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 (δ¹⁸O) in Greenland (Dansgaard et al., 1993; North Greenland Ice Core Project members, 2004) and also from several different proxies in Europe and in the Alps (e.g., Guiot et al., 1993; Genty et al., 2003; Spötl and Mangini, 2006). Two recent speleothems from a cave near Bern, Switzerland, also display these oscillations (Luetscher et al., 2015). These century- to millennia-long oscillations, called Dansgaard–Oeschger (DO) events, occurred repeatedly during Marine Isotope Stage 3 (MIS3, 60–30 ka BP) (Railsback et al., 2015). 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 (e.g., Keller and Krayss, 2005a; Ivy-Ochs, 2015; Beckenbach et al., 2014). This twofold maximum is also observed elsewhere in the Alps, on the Tagliamento glacier in the southern Alps in Italy (Monegato et al., 2007) and in Austria (van Husen, 1997). 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 Krayss, 2005b). 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 12 Comparison of geomorphic reconstruction with extent of temperate ice for simulation s5. (a) Geomorphic ice surface elevation of Benz-Meier (2003). 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 km², respectively, against 16 400 km²), 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 (Ellwanger et al., 2011). 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. Benz-Meier (2003) and Keller and Krayss (2005b) estimate an ice volume of ∼6500 km³, similar to simulation s1 (5500 km³) but significantly less than simulations s2–s5 (10 000 to 15 000 km³; see Table 2). Thicker ice for the Rhine glacier was also obtained by Becker et al. (2016) and Seguinot et al. (2018) 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}\left(\mathrm{100}{\mathrm{m}}^{-\mathrm{1}}\right){\mathrm{a}}^{-\mathrm{1}}$ and 0.025 $\mathrm{m}\mathrm{100}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$ for the ablation and accumulation zones, respectively. This value of ablation gradient is identical to the value obtained by Haeberli and Penz (1985) 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}\left(\mathrm{100}{\mathrm{m}}^{-\mathrm{1}}\right){\mathrm{a}}^{-\mathrm{1}}$ yields an ablation rate at the terminus (using 550 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ as a reference) of 0.65 m a⁻¹. 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⁻¹. The ablation rate can be converted back to a temperature with the model of Reeh (1991). Using a PDD factor of 6 mm PDD⁻¹ 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 (Reeh, 1991, Fig. 2). Higher PDD factors, as perhaps are found in Alpine glaciers (Hock, 2003), 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 Benz-Meier (2003) and Keller and Krayss (2005b) 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⁻¹, into 2.1 and 3.9 m a⁻¹ of melt at the terminus. Higher PDD factors would increase these values. Using Benz-Meier (2003) and Keller and Krayss (2005b) 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}\left(\mathrm{100}{\mathrm{m}}^{-\mathrm{1}}\right){\mathrm{a}}^{-\mathrm{1}}$ for Benz-Meier and Keller and Krayss, respectively, instead of 0.1 $\mathrm{m}\left(\mathrm{100}{\mathrm{m}}^{-\mathrm{1}}\right){\mathrm{a}}^{-\mathrm{1}}$ for simulation s1. For the accumulation area, using rough estimates of precipitation of Haeberli and Penz (1985) and the area of the accumulation zone based on hypsometric curves yields gradients of 0.034 and 0.06 $\mathrm{m}\left(\mathrm{100}{\mathrm{m}}^{-\mathrm{1}}\right){\mathrm{a}}^{-\mathrm{1}}$ for Benz-Meier and Keller and Krayss, respectively, instead of 0.025 $\mathrm{m}\left(\mathrm{100}{\mathrm{m}}^{-\mathrm{1}}\right){\mathrm{a}}^{-\mathrm{1}}$ for simulation s1. Gradients calculated on the basis of Benz-Meier (2003) or Keller and Krayss (2005b) 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 Machguth et al. (2016) 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}(\mathrm{100}\mathrm{m}{)}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$ in the very north and 0.45 m (100 m) a⁻¹ 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 (Fountain et al., 2006; Machguth and Cohen, unpublished) indicate values between 0.02 and 0.05 $\mathrm{m}(\mathrm{100}\mathrm{m}{)}^{-\mathrm{1}}{\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.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 (Fountain et al., 2006). 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 (Haeberli and Schlüchter, 1987), 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 Benz-Meier (2003), Keller and Krayss (2005b), or Bini et al. (2009). 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 (e.g., Benz-Meier, 2003; Kelly et al., 2004; Bini et al., 2009; Wirsig et al., 2016), does not represent the ice surface but an englacial cold–temperate ice transition layer that indicates the limit of glacial erosion (Florineth, 1998). 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 (e.g., Fabel et al., 2002; Briner et al., 2003; Ballantyne and Stone, 2015). The same situation is also believed to exist in some places under Antarctica (Creyts et al., 2014). 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 Benz-Meier (2003) 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 Benz-Meier (2003), 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.

Benz-Meier (2003)

Table 4 Annual 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 Benz-Meier (2003) 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⁻¹, 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 Benz-Meier (2003) 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⁻¹, about 30 to 80 times greater than computed for the geomorphic reconstruction using Stokes flow (<1 m a⁻¹; 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 Haeberli and Penz (1985) that included sliding yielded a value of 25 m a⁻¹.

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 (Beckenbach et al., 2014) 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 13 Three-dimensional flow lines for simulations (a) Benz-Meier (2003) 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 Benz-Meier (2003) (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 km³) and the highest rate (s5, 4.8 km³). 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 km³, 56–67 % depending on the simulation) than the Linth glacier (0.11 to 1.1 km³). Fluxes computed for the geomorphic reconstruction of Benz-Meier (2003) 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 Keller and Krayss (2005b) yielded a value of 1.86 km³, a value higher than the driest simulation (s1, 1.1 km³) but substantially smaller than the next two wetter simulations (s2 and s3, 3.3 and 3.1 km³, 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 14 Extent and thickness in meters (blue to red color map) of temperate basal ice for (a) Benz-Meier (2003) 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 (Blatter and Haeberli, 1984; Haeberli and Penz, 1985), 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 (Pohjola and Hedfors, 2003; Lüthi et al., 2015) and in conceptual models (Krabbendam, 2016). 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 15 Vertically integrated strain heating in watts per square meter for (a) Benz-Meier (2003) 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⁻² for simulations s2–s5. In comparison, geothermal heat flux is of the order of 0.1 W m⁻² and does not play a significant role in determining basal temperature. Geothermal heat flux, which varies between 0.06 and 0.12 W m⁻² 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 (Alley, 1993), 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 16 Effect 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⁻² 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⁻¹ with even slower speeds in the lobe (average surface speed of 25 to 80 m a⁻¹). 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.