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 km2. As part of an integrative study about the safety of repositories for radioactive waste under ice age conditions in Switzerland, we modeled the Rhine glacier using a fully-coupled, three-dimensional, transient, thermo-mechanical Stokes flow 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 (TJuly 0 °C at the glacier terminus) and dry ( 10 to 20 % of today’s precipitation) climatic conditions. Forcing the numerical simulations with warmer and wetter conditions that better matched LGM climate proxy records yielded a glacier on average 500 to 700 m thicker than geomorphic reconstructions. Mass balance gradients also controlled ice velocities, fluxes, and sliding speeds. These gradients, however, had only a small effect on basal conditions. All simulations indicated that basal ice reached the pressure melting point over much of the Rhine and Linth piedmont lobes, and also in the glacial valleys that fed these lobes. Only the outer margin of the lobes, bedrock highs beneath the lobes, and Alpine valleys at high elevations in the accumulation zone remained cold based. The Rhine glacier was thus polythermal. Sliding speed estimated with a linear sliding rule ranged from 20 to 100 m a−1 in the lobes, and 50 to 250 m a−1 in Alpine valleys. Velocity ratios (sliding to surface speeds) were > 80 % (lobes) and 60 % (valleys). Basal shear stress was very low in the lobes (0.03–0.1 MPa), 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 W m−2 but on average 0.03 to 0.2 W m−2. Basal friction acted as a heat source at the bed of about 0.02 W m−2, 4 to 6 times less than the geothermal heat flow which is locally high (up to 0.12 W m−2). In the lobes, despite low surface slopes and low basal shear stresses, sliding dictated main fluxes of ice which closely followed bedrock topography: ice was channeled in between bedrock highs along troughs, some of which coincided with glacially eroded overdeepenings. These sliding conditions may have favored glacial erosion by abrasion and quarrying. Our results confirmed general earlier findings but provided more insights into the detailed flow and basal conditions of the Rhine glacier at the LGM. Our model results suggested that the trimline could have been buried by a significant thickness of cold ice. These findings have significant implications for interpreting trimlines in the Alps and for our understanding of ice-climate interactions. Abstract. At the Last Glacial Maximum (LGM), the Rhine glacier in the Swiss Alps covered an area of about 16,000 km 2 . 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 fully-coupled, three-dimensional, transient, thermo-mechanical Stokes ﬂow 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 July ∼ 0 ◦ C at the glacier terminus) and dry ( ∼ 10 5 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, ﬂuxes, 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 10 beneath the lobes, and Alpine valleys at high elevations in the accumulation zone remained cold based. The Rhine glacier was thus polythermal. Sliding speed estimated with a linear sliding rule ranged from 20 to 100 m a − 1 in the lobes, and 50 to 250 m a − 1 in Alpine valleys. Velocity ratios (sliding m − 2 Basal , to the locally to ). the basal was in between bedrock along may have Our results conﬁrmed general earlier ﬁndings but provided more insights into the detailed ﬂow and basal conditions of the Rhine glacier at the LGM. Our model results suggested that the trimline could have been buried by a signiﬁcant thickness of cold ice. These ﬁndings have signiﬁcant implications for interpreting trimlines in the Alps and for our understanding of ice-climate interactions. ice ﬂow 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 20 gradients, and compared our simulations of ice thickness and ﬂow (where possible) to geomorphic reconstructions and earlier glaciological studies. Our numerical simulations conﬁrmed results of earlier studies but provided more details of the complex three-dimensional ice ﬂow 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 25 150 to 530 m a − 1 with even slower speeds in the lobe (average surface speed of 25 to 80 m a − 1 ). Sliding patterns mimicked surface speeds with 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 ﬂux and yielded complex patterns of ice ﬂow in the lobes. Basal shear stress in the lobe was smaller than 0.1 MPa for the high mass balance gradient simulations, and down to 0.025 MPa when prescribing low mass balance gradients. The latter value is similar to earlier calculations. Higher 30 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 corresponds to old interpreted as a cold-temperate ice surface above which glacial erosion did not occur, leaving no erosional features on mountain valley ﬂanks. Our numerical results indicate that transection glaciers such as the Rhine glacier have complex three-dimensional ﬂow patterns that are difﬁcult to reduce to two dimensions. The coupling of ice ﬂow 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 precision needed to fully understand 5 the coupling between ice ﬂow, climate, and basal thermal conditions. Future models that include the effects of permafrost on basal sliding, and how ice loading and distribution inﬂuence patterns of ground water ﬂow are needed to better understand conditions at the LGM and the characteristics of the subsurface during ice-age conditions.


Introduction
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 5 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;Keller, 1988;Beckenbach et al., 2014). Erratics, tills and other sediment 10 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 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 25 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 are physically not consistent when sliding is an important component of flow, neglected important stresses (longitudinal and transverse) relevant to the complex ice flow patterns of transection glaciers 30 (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 35 3 The Cryosphere Discuss., https://doi.org /10.5194/tc-2017-204 Manuscript under review for journal The Cryosphere Discussion started: 12 October 2017 c Author(s) 2017. CC BY 4.0 License.  (Pietsch and Jordan, 2014). Outline of maximum glacial extent at LGM (Bini et al., 2009) is shown in blue. 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 5 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 estimate more precisely the flow of the Rhine glacier at the LGM. 10 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 mid-and high-latitude regions can be affected by impacts from future ice-age conditions. Several countries have developed programs to investigate 15 potential future ice-related environmental changes and their effects at depth (Fischer et al., 2015). The main concerns are deep 4 The Cryosphere Discuss., https://doi.org /10.5194/tc-2017-204 Manuscript under review for journal The Cryosphere Discussion started: 12 October 2017 c Author(s) 2017. CC BY 4.0 License. erosion by glaciers or ice sheets, penetration of permafrost to great depth, 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 5 of repositories for radioactive waste in northern Switzerland (see Fig. 3).
Here, we use a fully three-dimensional, numerical, thermo-mechanical 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 solu- 10 tions, 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 15 speeds, basal temperatures and shear stresses from which ice fluxes, flow patterns, and areas of temperate basal condition can be derived.

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, 20 open source code Elmer/Ice (Gagliardini et al., 2013). Despite its significant computational cost, Stokes flow is more accurate than shallow ice/shallow shelf approximation in zones where gradients in stresses are significant like in complex transection glaciers found in Alpine settings and were 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 25 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 south-east 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 30 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.

5
The Cryosphere Discuss., https://doi.org /10.5194/tc-2017-204 Manuscript under review for journal The Cryosphere Discussion started: 12 October 2017 c Author(s) 2017. CC BY 4.0 License. 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.
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 Paul and Linsbauer, 2012), 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. For the 5 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 (Jäckli, 1962;Jäckli, 1970) and also of Keller and Krayss (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 10 and the ice surface.

Thermomechanical model
Stokes equations describe the mass and momentum balance for a viscous fluid. Together with the equation for energy, this system forms a thermo-mechanical 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.

Stokes equations
The conservation of mass for an incompressible fluid (divergence-free velocity field) is where ∇· indicates the divergence. The conservation of momentum is 10 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 where d = ( 1 2 tr(D 2 ) is the effective strain rate, d o is the critical strain rate, D = 1 2 (∇v) + (∇v) T 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 15 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.

The heat equation 20
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 Geothermal heat flow in the area of the LGM Rhine glacier is highly variable ranging from 0.06 to 0.12 W m −2 with low values in the high Alps and high values in geothermally active regions of the Swiss Plateau. We use the heat flow data of Medici and Rybach (1995) obtained from numerous measurements of temperature gradients in boreholes and of rock thermal 5 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, non-linear way. Correcting for these effects is difficult and is not done in Medici and Rybach (1995). Neither is it done here. It would require, at a minimum, to impose 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.

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). Sub-freezing sliding can also occur (e.g., Shreve, 1984) but is small and negligible for ice fluxes or erosion for the time scale considered here (100s to 1000s years). It is neglected here. A commonly-used approach for the sliding speed is to link it to the basal shear stress, i.e., where v s is the sliding speed vector, τ b the basal shear stress vector, m is an exponent, and C is 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 5 we assume m = 1 (linear sliding).
Ice sliding over the substrate occurs when temperatures reach the melting temperature. Below it 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), where C o and C 1 are the sliding coefficients for temperate and cold conditions, respectively, and T = T − T pmp where T pmp is 10 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.

Thermal boundary condition
The boundary condition for the basal temperature is (Seddik et al., 2012):

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

20
Solving for ice velocity, temperature, and the elevation of the ice surface requires three types of boundary condition 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, (iii) and the ice temperature.

Stress-free surface
The glacier surface is a free surface in contact with the atmosphere. This surface supports no shear stress. This is simply 25 expressed as t · n = 0, where t is the stress vector at the ice surface and n is the unit normal pointing outward. 9 The Cryosphere Discuss., https://doi.org /10.5194/tc-2017-204 Manuscript under review for journal The Cryosphere Discussion started: 12 October 2017 c Author(s) 2017. CC BY 4.0 License.

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 where ζ is the glacier surface elevation andḃ 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/ablation function.

Surface mass balance
For the accumulation/ablation function (ḃ in Eq. (9)), we choose a simple parameterization represented by two mass balance gradients, one for the accumulation area, one for the ablation area, and a maximum threshold value for the maximum accumu-

10
where z is the elevation of the ice surface , β acc and β abl are the accumulation and ablation mass balance gradients, respectively, andḃ max acc 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 (next section), we have effectively decoupled the mass balance from the energy. A rate of ablation based on the number of positive degree-day (PDD model) would have been 15 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;Jouvet et al., 2017). Furthermore, PDD factors needed to compute surface melt rates are known to vary substantially (cf. Braithwaite, 1995;Hock, 2003;van den Broeke et al., 2010) and choosing suitable factors for LGM conditions is challending. The same is true for accumulation which requires knowing patterns of rainfall and temperature. Since LGM 20 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ḃ max were selected to represent a range of dry to wet climates while values of β abl were chosen to be either in agreement with earlier estimates for the LGM Rhine 25 glacier (Haeberli and Penz, 1985) or with present-day Greenland values (Machguth et al., 2016) measured in regions where temperatures and precipitation are similar to estimated LGM conditions in northern Switzerland.

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 done in many ice modelling 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 reconstruction 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, 1983Haeberli, , 2010. We assume 5 that changes in surface temperature with elevation are linearly related to a lapse rate, γ a , where T ela and z ela are the temperature and elevation at the equilibrium line, respectively. In all model simulations γ a =  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).

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 between 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 25 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 = −ḃ z/H is, for the accumulation zone (ḃ > 0), where [dT /dz] bed is the gradient of temperature at the bed and z 2 = 2α T H/ḃ. For ablation (ḃ < 0), the temperature distribution where D() 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 10 computed temperature sometimes exceeds the melting temperature when the ice is thick. This is not surprising since nowhere in the model of Robin (1955) there is information about a maximum temperature or about a phase transition. Limiting the temperature to the melting temperature is done 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 geother-15 mal 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.

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). Scenario s1 has the smallest accumulation 20 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. Other parameters (basal topography, glaciological parameters) 25 are identical for all simulations. All input parameters and key computed quantities for these simulations are summarized in Table 2. For comparison, the table also includes data for two geomorphic reconstructions: Benz-Meier (2003) and Keller and Krayss (2005b).

12
The Cryosphere Discuss., https://doi.org /10.5194/tc-2017-204 Manuscript under review for journal The Cryosphere Discussion started: 12 October 2017 c Author(s) 2017. CC BY 4.0 License. Table 2. Summary of numerical (this study) and selected geomorphic reconstructions. zela is elevation of equilibrium line. βabl and βacc are the ablation and accumulation mass balance gradients, respectively.ḃmax is the maximum accumulation rate. bacc is the average accumulation rate. bnet is the glacier net mass balance.ḃterm is the ablation rate at the terminus. AAR is the accumulation area ratio. A is the glacier area, V its volume. ts is the simulated time. Numerical simulations s1 through s5 are ranked from driest to wettest based on the value of the average accumulation rate, bacc. Initial conditions for the ice surface are: the geomorphic reconstruction of Benz-Meier (2003) for s1, a simulation with βacc = 0.1 m (100m) −1 a −1 , βabl = 0.2 m (100m) −1 a −1 , and zELA = 1200 m that ran for 440 years for s2, and a simulation with and zELA = 900 m that ran for 907 years for s3, s4, and s5.
Reconstruction 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 (ḃ max ). Model outputs are the average net accumulation (b acc ), the average net glacier balance (b net ), the specific balance rate at the terminus (ḃ 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 5 and volume obtained from field mapping and inferences of ice surface elevation contours, AAR,ḃ term , andb net assumed to be zero (the glacier is in a state of dynamic equilibrium).
For the geomorphic reconstructions, the termḃ 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 positive degree days (PDD) using the model of Reeh (1991), and then 10 by multiplying the PDD by a PDD factor, here taken equal to 6 mm PDD −1 (Braithwaite, 1995). 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 5 calculated for Keller and Krayss (2005b) Table 2) is higher than the value cited in Keller and Krayss (2005b) (0.30 m a −1 ) obtained from numbers cited in Haeberli (1991) that assumed LGM precipitation was equal to 20% of today's value (1.5 m a −1 ).

Numerical solutions
The open source software Elmer/ice (Gagliardini et al., 2013) is used to solve the set of equations and their boundary conditions 10 using the finite element method. The general 3-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 15 a steady state. Simulations were stopped due to computational time constraints. 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ḃ net = 0.02 m a −1 , see Table 2) and simulation s3 seems close to it (ḃ net = 0.19 m a −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 (e.g., Levermann and Winkelmann, 2016). Whether this applies to our 25 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).
14 The Cryosphere Discuss., https://doi.org /10.5194/tc-2017-204 Manuscript under review for journal The Cryosphere Discussion started: 12 October 2017 c Author(s) 2017. CC BY 4.0 License.  (2003) and (b-f) simulations s1 through s5 (Table 2). Black contours are every 500 meters. The 2500 m contour is shown in yellow. The magenta contour indicates the equilibrium line altitude. 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.
Of these two simulations, simulation s3 is closer to the LGM margin with an ice surface area of 14,553 km 2 while simulation 30 s1 is only 12,389 km 2 (the geomorphic reconstructions yielded values around 16,000 km 2 ). The ice surface area of simulation s5 is nearly identical to the geomorphic reconstruction (16,400 km 2 ). Simulation s4 is slightly larger (17,706 km 2 ); s2 slightly smaller (15,813 km 2 ). Their configurations, however, are significantly different with most Alpine peaks and the Hörnli ridge ice covered, 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 5 move north, beyond the LGM margin. 15 The Cryosphere Discuss., https://doi.org /10.5194/tc-2017-204 Manuscript under review for journal The Cryosphere Discussion started: 12 October 2017 c Author(s) 2017. CC BY 4.0 License.
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)

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 5 numerical simulations and for the Benz-Meier (2003) and Keller and Krayss (2005b) geomorphic reconstructions.

15
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 AAR can be read directly from the graph, the AAR being the intersection of the ELA with the hypsometric curve. All numerical reconstructions have 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 20 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.

Basal conditions
Basal temperature and sliding speed are two key variables that describe basal conditions. Basal temperature is shown in Fig. 8 as the difference between modelled ice temperature and the pressure melting point calculated as a function of ice pressure at the bed. The small range shown in Fig. 8 (Table 2). Yellow contour indicate the extent of ice that is within 0.05 • C of the pressure melting point.
simulation based on the geomorphic ice surface reconstruction (Fig. 8a), the bed is assumed temperate as in Benz-Meier (2003).

5
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).

5
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 tem- 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.
Following patterns of surface speed (Fig. 7), sliding speed is highest in the Alpenrhein in simulation s5 (Fig. 9f) reaching 5 about 200 m a −1 . Maximum sliding speed diminishes for simulations which have smaller ice throughflow and surface slopes.
In simulation s1 which has the least amount of ice, maximum sliding speed is less than 100 m a −1 at the Sargans diffluence. In the lobes, sliding speed generally diminishes quickly as ice radially spreads into the lowlands but sliding speed patterns there also follow topographic patterns of troughs and bedrock highs. For example, sliding speed is relatively high (150 m a −1 ) in the two main branches (Limmat and Glatt valleys) of the Linth lobe. In the larger Rhine lobe, sliding is significant along the Lake Constance trough and the Thur valley near Frauenfeld (see Fig. 3 for locations), with speeds in excess of 100 m a −1 .
These patterns are more pronounced for simulation s5 which has the highest ice flux.
Basal shear stress, another indicator of basal condition, is shown in Fig. 10 and is calculated using Eq. 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 (Braedstrup et al., 2016), and with 15 local measurements beneath glaciers (Cohen et al., 2000;Iverson et al., 2003;Cohen et al., 2005). 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.
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 20 greater (80% of surface speed) in the lobes than in the main trunks of valleys (50-70%). In the lobe, ice flows almost like a plug with little 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 where 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 25 main Alpine valleys, sliding is about 50 to 60% 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.

Steady state
Numerical solutions simulated at least 1500 years of glacier evolution and more than 3000 years in one case (see Table 2).  (Table 2). 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).
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 condition at the LGM. The climate record prior to the LGM indicates large oscillations as shown by the oxygen isotopic record (δ 18 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., 5 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 timescale oscillations, called Dansgaard-Oscheger (DO) events, lasted through the Marine Isotope Stage 3 (MIS3, 60-30 ka BP) until the LGM (Railsback et al., 2015). The climate 23 The Cryosphere Discuss., https://doi.org /10.5194/tc-2017-204 Manuscript under review for journal The Cryosphere Discussion started: 12 October 2017 c Author(s) 2017. CC BY 4.0 License.   Keller and Krayss, 2005a;Ivy-Ochs, 2015;Beckenbach et al., 2014). This two-fold 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, 5 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

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 2 , respectively, against 16,400 km 2 ), 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 5 reconstruction. The former is due to overall thicker ice for these two simulations. The latter may be due to using the presentday 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 10 lobe.
In addition to glacial extent, glacial volume may help discriminate between numerical results. Except for simulation s1, all other simulations indicate a thicker glacier (by at least 500 meters) than the geomorphic reconstructions. Benz-Meier (2003) and Keller and Krayss (2005b) estimate an ice volume of ∼ 6500 km 3 , similar to simulation s1 (5500 km 3 ) but significantly less than simulations s2-s5 (10,000 to 15,000 km 3 , see Table 2). Thicker ice for the Rhine glacier was also obtained by Becker 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 m (100 m −1 ) a −1 and 0.025 m 100 m −1 a −1 for the ablation and accumulation zones, respectively. This value of ablation 20 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. These mass balance gradients yield an ablation rate at the terminus of 0.65 m a −1 and an average accumulation rate of 0.19 m a −1 , numbers representative of an extremely cold and dry climate. Converting these numbers to temperature and precipitation using PDD (Reeh, 1991) yields a LGM July temperature of about 0 • C and a LGM precipitation between 10 25 and 20% of today.
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 of 3.2 and 7 • C, respectively, which translates, using the positive degree day model of Reeh (1991) with a PDD factor of 6 mm PDD −1 , into 2.1 and 3.9 m a −1 of melt at the terminus. Using their values of ELA (see Table 2), these numbers yield ablation gradi- 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, simulation that best match the Rhine glacier reconstructed from geomorphic interpretation, is smaller than any present-day values measured in Greenland. Furthermore, the gradients of simulation  (Fountain et al., 2006;Machguth and Cohen, unpublished) indicate values between 0.02 and 0.05 m (100 m) −1 a −1 , the latter is only one-half of the value used in simulation s1. Estimated melt 10 rates at the terminus of the Greenland Ice sheet and of local glaciers at an altitude similar to the margin of the Rhine glacier lobe range between 1 and 5 m a −1 , 1.5 to 7 times larger than the value used in simulation s1. The marginal melt rates 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 15 (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 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, 20 up to 700 m more, than the geomorphic map of Benz-Meier (2003), Keller andKrayss (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., 25 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, 30 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 26 The Cryosphere Discuss., https://doi.org /10.5194/tc-2017-204 Manuscript under review for journal The Cryosphere Discussion started: 12 October 2017 c Author(s) 2017. CC BY 4.0 License. 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 5 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 near 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 10 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 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 m a −1 to 535 m a −1 , values that reflect overall small mass balance gradients. Highest flow speeds are obtained for simulation s5 which has the highest rate of accumulation. These surface speeds are significantly smaller than the maximum surface speed estimated from the geomorphic reconstruction. Differences between the geomorphic (Fig. 7a) and numerical reconstructions (Fig. 7b-f) are due to the evolution of the ice free surface for (c) Same as (a) with temperate-cold transition of (b) superimposed. All elevations are in meters.
using Stokes flow (< 1 m a −1 , see V 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 −1 . Figure 7 shows surface velocity but also included ice flowlines 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 5 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).

Flow patterns and fluxes
Results indicate that the Rhine lobe is made up of ice originating from the Valserrhein, Hinterrhein, Landquart, and Ill Additional insight about flow patterns can be gained by generating tracers randomly in a basin and following their paths downstream. Figure 13 shows ice flowlines colored according to the source region of the ice. Flowlines are calculated from 5 the three-dimensional velocity field. For the geomorphic reconstruction (Fig 13a), Ice thickness also influences flow paths with thicker ice allowing movement of ice across high passes and ridges in the Alps.

20
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. 25 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 3 ) and the highest rate (s5, 4.8 km 3 ). The Rhine glacier, because of its large basin area, brings most ice to the Rhine lobe (about 80%), with little variation between the simulations (s1-s5), the remaining ice in the Rhine lobe coming from the Ill glacier. For the Linth lobe, The Walensee glacier 30 brings slightly more ice (0.27 to 1.4 km 3 , 56-67% depending on the simulation) than the Linth glacier (0.11 to 1.1 km 3 ).
Fluxes computed for the geomorphic reconstruction of Benz-Meier (2003) Table 4). Background color is elevation of basal surface.
glaciers in the Linth lobe. Given the vastly different basin sizes of these glaciers, these numbers do not appear to be realistic. 35 An estimate of ice flux entering the Rhine lobe by Keller and Krayss (2005b) yielded a value of 1.86 km 3 , a value higher than the driest simulation (s1, 1.1 km 3 ) but substantially smaller than the next two wetter simulations (s2 and s3, 3.3 and 3.1 km 3 , respectively).

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 where 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.

5
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  (Table 2). Background green to brown color scheme shows basal topography.

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 5 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,   (Table 2).
probably because of flow constriction associated with the valley width and 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 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 there where 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.

5
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, bedrock ridges and bumps, confluences, are zones where 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 10 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 to 0.1 MPa for thin and thick-ice simulations, respectively.
In one simulation we turned on basal friction (see Eq. (8)). Figure 16 shows the effects of basal friction on basal conditions. Basal friction generates additional heat at the glacier bed (about 0.02 W m −2 on average), increasing basal temperature to the melting temperature everywhere in the lobe, increasing the sliding speed by a few 10's of m a −1 everywhere on the Swiss

15
Plateau. This 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.

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 uneroded gravel deposits, several hundreds to one 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. Simulations that used a less extreme cold-dry climate yielded an ice thickness 500 to 700 meters more than the interpreted as a cold-temperate ice surface above which glacial erosion did not occur, leaving no erosional features on mountain valley flanks.
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 precision needed to fully understand 5 the coupling between 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 conditions at the LGM and the characteristics of the subsurface during ice-age conditions. grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the program Investissements d'Avenir supervised by the French Agence Nationale pour la Recherche. Thanks are due to Olivier Gagliardini, Thomas Zwinger, and Martina Schäfer for help with Elmer/Ice.