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

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

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

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,

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

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

Table 1Model parameter values.

Download Print Version | Download XLSX

2.4 Basal boundary conditions

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.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}$),

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

where 𝒟() is the Dawson integral

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 2Summary 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.

Download Print Version | Download XLSX

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.1 Ice extent, surface elevation, and ice thickness

Figure 4Ice 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.

Download

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 5Ice thickness in meters for (a) Benz-Meier (2003) and (b–f) simulations s1–s5.

Download

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 6Glacier hypsometric curves for geomorphic and numerical reconstructions of the Rhine glacier system. Vertical colored lines show the value of the AAR on the x axis for each reconstruction.

Download

3.3 Ice surface speed

Figure 7Ice 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.

Download

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.

Table 3Summary 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.

Download Print Version | Download XLSX

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 8Basal 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.

Download

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 9Ice sliding speed in meters per year for (a) Benz-Meier (2003) (no sliding assumed) and (b–f) simulations s1–s5.

Download

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 10Basal shear stress in megapascals for (a) Benz-Meier (2003) (steady-state Stokes solution) and (b–f) simulations s1–s5.

Download

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

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

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

Download

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 12Comparison 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.

Download

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.

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

Download Print Version | Download XLSX

4.3 Rhine glacier dynamics

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 13Three-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.

Download

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 14Extent 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.

Download

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 15Vertically 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.

Download

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

Download

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.