A study of the Würm glaciation focused on the Valais region (Alps)

Abstract

During the Last Glacial Maximum (LGM), the glaciation in the European Alps reached maximum ice extent. We already simulated the steady states of the Alpine ice coverage for several climate drivers in Becker et al. (2016) and heighten in this article such studies for the Swiss Valais region. To this end, we employ the Parallel Ice Sheet Model (PISM), which combines the shallow ice approximation (SIA) with basal sliding elements of the shallow shelf approximation (SSA), and subject this model to various external driving mechanisms. We further test the sensitivity of this kind of the ice coverage in the Valais region to a temporally constant climate and to monotonic ice sheet build-up from inception to steady state as well as to the Dye 3 temperature driving during the past 120 000 years. We also test differences in the precipitation patterns exerted to the northern and southern catchment areas of the Rhone and Toce rivers to possible transfluence changes in ice from the northern to the southern catchment areas and vice versa. Moreover, we study the effect of the ice deformability and estimate the removal up to 1000 m of sediment in the Rhone Valley and study the removal of rock hindering the flow through the valley cross section at the knee of Martigny. All these studies took place because of a discrepancy in the ice height prediction of the modelled ice sheet with its geomorphologically reconstructed counterpart with proxy data obtained by Bini et al. (2009) as well as a difference in ice height between the two of up to 800 m. Unfortunately, all the scenarios in the model do not sufficiently reduce this discrepancy in the height prediction and the geomorphological reconstruction. The model results have discovered an ice dynamical discrepancy with the land map in Bini et al. (2009).

Introduction

Preamble

This paper is an excerpt of a dissertation (in the German language) (Becker, 2017) that focuses on the numerical modelling of the glaciation of the Alps during the Last Glacial Maximum. The dissertation first concentrates on the more local response of the Valais region in the Swiss Alps to refined external climate forcing. The current paper is devoted to this behaviour in the glaciation of the Valais region but goes slightly beyond the interpretation in the second part of Becker (2017).

In the ensuing discussion of the results, particular locations in the Valais region are defined by nomenclatures of today's custom, even when their locations may have been covered by ice, i.e. did not exist in earlier times. As an example, Brig means the location on the land map which today is occupied by the town Brig.

Historical background

The traces left behind by the past glaciations of the Alps contribute in an essential way to the current understanding of the glacial extent and the dynamics of the Alpine ice sheet as well as to the implicit comprehension of the climate variations during the past four ice ages of the European Alpine region. The reconstruction of the glaciation in this region is in our opinion particularly apt when being applied to the maximum glaciation extent of this region, because the left-behind traces are often optimally recognized, since these traces are less overshadowed by variations in advances and retreats of glaciers than otherwise. As noted in Becker et al. (2016), subsequently simply referred to as (I), terminal moraines are particularly apt to infer information on the lateral extent of the glaciated area. Similarly, glacial erratics (in this article erratic blocks or boulders will be called “erratics” for brevity.) function likewise; they may additionally provide information on the flow dynamics of the glacial sheets. Conversely, in order to determine the past vertical glacial extent, boundaries of glacier polish seem to be the best markers with which to identify the past glacial extent.

With the aid of this evidence, a better understanding of the last glaciation of the Alpine region has already been developed over past decades. The inferred understanding, gained by past geomorphological studies of the glaciation of the Alps, described the latter as a glacial-flow net dictated by the relief of the underground, similarly to today's glaciation of southern Alaska. Jäckli (1970) systematically evaluated and interpreted the glacial evidence. This led to the first map of the LGM glaciation of Switzerland. However, the interpretation of the glaciation of Switzerland shown in this map was rather strongly influenced by the glacial-flow net idea.

Later, in Bini et al. (2009), a revised version of the Jäckli map was published. In this revision, as well as in articles by Florineth (1998), Florineth and Schlüchter (1998), (2000) and Kelly et al. (2004), the glacial-flow-net-like inferred topography of the free ice surface was questioned. Moreover, indications were presented which hinted at a dome-like glacial surface in the regions of the present Rhone Glacier (Rhone ice dome) and in the Surselva (Anterior Rhine dome) as well as in the region of the Bernina mountain (ice dome over Upper Engadin); see e.g. Jäckli (1970), Bini et al. (2009) or Florineth (1998). Characteristically, in these articles the glaciation of the western Alps is still based on the idea of a glacial-flow net, which, apart from the Rhone ice dome, is particularly nourished by the ice fields that are embedded at the valley ends. Apart from the Rhone ice dome of the Aletsch region, the ice field in the Mattertal and the Mont Blanc region are identified here as the principal accumulation centres. The results of the models in (I), however, point to an ice-dome-like free-surface topography of the glaciated area of the western Alps, particularly in the Upper Valais. This interpretation holds, in particular, for conditions for which the lateral extent of the ice at the LGM is to be reproduced.

The morphology of the glaciation is chiefly determined by the ice thickness. It is indeed the ice thickness that determines whether the peaks are covered by ice and so enables ice caps or ice domes to be formed, from which the surficial ice transport occurs essentially approximately radially or down lines of steepest surface descent, or whether the ice height lies below the mountain peak, and the ice flow is guided by the mountain flanks.

In (I), it was shown that the geomorphologically reconstructed lateral extent of the LGM glaciation can be modelled by the available glacio-mathematical descriptions that were developed in the last 40–50 years. Such computations deduced from fundamental classical physical laws (balances of mass, momentum and energy) are subject to simplifications in order to reduce actual CPU times in realistic modelling scenarios (for specifics see (I) or Hutter, 1983; Greve and Blatter, 2009; Hutter et al., 1986; Hutter and Wang, 2016; the PISM authors, 2015 and others). The result is, however, that even in situations where the lateral extent of a computed LGM ice mass is in satisfactory agreement with observed proxy data, the vertical extent, i.e. the glacial height of the computed ice dome, deviated considerably from information obtained with proxy data. Particularly strong deviations from the geomorphologic reconstruction are observed in high-altitude inner-Alpine valleys such as the Surselva and the Rhone Valley in the Valais region. As the modelled free ice surface in the Surselva lies about 500 m above that of the geomorphologically reconstructed geographical map of Bini et al. (2009), the corresponding deviation in the Upper Rhone Valley is as large as 800 m.

As a further novelty compared to the Jäckli map, the map of Bini et al. (2009), which is based on the improved study of Kelly et al. (2004), evidences a transfluence of the ice from the Upper Rhone Valley over the Simplon Pass into Val Diverdo and thus into the catchment area of the river Toce.

In this paper, it is our intention to search for the vertical extent of the ice coverage of the Upper Rhone Valley and the surrounding mountain ranges at LGM; this roughly corresponds to the area of today's glaciers in Canton Valais. This will be our focus in this article, since the deviation of the modelled ice height (of the free surface) from the geomorphological ice height in previous studies has been shown to be particularly large. In the frame of this context we shall also discuss the morphology of the free surface at the LGM of the glaciation and the configuration of the glacier flow system. In particular, our focus will also concentrate on the region of the Simplon Pass with special attention on the geomorphological evidence of a possible north–south transfluence.

Modelled region

The modelled area will be restricted to the region of the Upper Valais. In order to reduce possible unrealistic perturbations due to boundary effects of the model, a distinctly larger numerical region will be chosen as our basic computational region. This region is between today's Lake Geneva in the west (E6${}^{\circ }$40${}^{\prime }$), Surselva in the east (E8${}^{\circ }$50${}^{\prime }$), the region of the Napf in the north (N47${}^{\circ }$) and the southern boundary of the Alps (N45${}^{\circ }$40${}^{\prime }$). Of primary interest is the Upper Rhone Valley with the surrounding mountain ranges of the Bernese and Valais Alps. Both sides of the Upper Rhone Valley (Goms) are dominated by high mountain chains, the peaks of which surpass the 4000 m a.s.l. mark. In the north, the Finsteraarhorn (4274 m a.s.l.) and the Aletschhorn (4193 m a.s.l.) form the highest peaks. At their feet lies the Grosser Aletsch glacier, by area the largest glacier of the Alps, and quite a number of additional ice streams. The Alps of the Valais extend to the south beyond the Rhone Valley; they comprise the largest number of peaks higher than 4000 m a.s.l. in the Alps and are dominated by the Monte Rosa Massif (4634 m a.s.l.).

In the ensuing study, particular attention will be focused on the region around the Monte Leone–Blinnenhorn chain (see Fig. 1, dashed black rectangular region). Since this rectangle is embedded into the much larger glaciated region, here the effects due to the boundaries are likely negligibly small. The Monte Leone–Blinnenhorn chain separates the Goms and the Binn Valley, Upper Valais, a catchment region of the river Rhone, from the southern valleys Val Formazza, Val Antigorio, Val Divedro and Val d'Ossola, a catchment regions of the river Toce (Tosa) in the Piedmont. At the western boundary of the Monte Leone–Blinnenhorn chain, between the Weissmies group (4017 m a.s.l.) and the Monte Leone (3553 m a.s.l.) lies the considerably less high Simplon Pass, which connects the Upper Rhone Valley in the surroundings of Brig and following the Saltina river to the top of the pass with the Gondo Valley, which merges into the Val Divedro. The Val Divedro discharges into the Val d'Ossola, as do the transitional stretches of the Monte Leone–Blinnenhorn chain. Contrary to these, the Simplon Pass is easily accessible by road for field studies.

Extent of the computational region. It consists principally of the region of the Rhone Valley between the Furka Pass and the beginning of Lake Geneva as ablation area of the ice age glaciation and the neighbouring mountain flanks as region of accumulation. The extent of the ice during the LGM after Bini et al. (2009) is shown with a blue background. The red dashed line shows the boundary of the weather as given and quantified in (I). Because at the boundaries of the simulated rectangular region, computational errors can arise due to reflections or glacier flow from outside into the rectangular region, our principal focus will be on the interior region of the Upper (Swiss) Rhone Valley. This region is shown once more in detail and comprises the region within the dashed black lines. Reproduced with permission from swisstopo (JA100120).

[Figure omitted. See PDF]

Methodology

For the numerical computations we employ the Parallel Ice Sheet Model (PISM) (see the PISM authors, 2015) in order to simulate the glaciation of the LGM in the discretized region of Sect. 2 (see also Fig. 2). In doing so the principal approach of the model strategy will be maintained. This means that the hybrid model described in (I) will again be employed. More specifically, a superposition of the shallow ice approximation (SIA) (Hutter, 1983; Hutter et al., 1986; Greve and Blatter, 2009) and the shallow shelf approximation (SSA) (Morland, 1987; Weis et al., 1999) will be employed, the latter of which will serve as a mechanism for a pseudo-plastic sliding law (Mohr–Coulomb friction).

Zoom of the region of the black-dashed rectangular region in Fig. 1, forming the physically relevant region around the Monte Leone–Blinnenhorn chain. This mountain chain today separates the two catchment areas of the rivers Rhone and Toce. Reproduced with permission from swisstopo (JA100120).

[Figure omitted. See PDF]

One particular proviso should be kept in mind when working with the PISM that applies SIA$+$SSA. The mathematical property of this is that the horizontal flow component of the ice is in any material point of the ice in the direction of the horizontal projection of the direction of steepest descent of the free ice surface. This means that in any real or hypothetical vertical ice borehole, the modelled horizontal component of the ice flow is in the direction dictated by the free-surface speed, irrespective of whether the considered ice point is at the base at which the basal sliding speed shares this direction. This property also holds in PISM if the basal sliding law is derived from the shallow shelf approximation (see e.g. Hutter, 1983; Weis et al., 1999). Therefore, the flow direction along a flow line holds for all depths of a given flow line. This indicates a possible limitation of the computed PISM results if observations point to a deviation of the real flow from this parallelism. In the latter case, at least the Stokes approximation (see Hutter and Wang, 2016) should be applied. However, the CPU times for computations with the Stokes model are far too large these days to be used in the present situation.

For the reference parameterization of the model used so far, starting from a regular grid in the initial position, the flows at the free ice surface were determined by using the lateral ($x,y)$-flow speed at the free ice surface and by supposing steady-state conditions. They are determined as the orthogonal trajectories of the level lines of the free surface. The flow lines of the momentary situation (at the chosen time slice) graphically emphasize the flow field of the ice at the chosen instant. In order to avoid an erroneous interpretation, it is worth noting that during the ice build-up these do not represent the ice particle trajectories. In particular, during a rapid ice build-up process the flow lines (streamlines in fluid mechanics) and particle trajectories are not the same lines. However, if the free ice surface changes slowly or under steady-state conditions, the two are close to each other or identical (for steady state); see e.g. Hutter and Wang (2016, p. 64–73). Because the streamlines guarantee a better representation, and the mentioned differences between the streamlines and particle trajectories are small and often negligible, we prefer the streamline representation.

The geothermal heat flow will be set equal to 75 W m${}^{-\mathrm{2}}$ This is the mean value suggested by Shapiro and Ritzwoller (2004). Moreover, the angle of internal friction of the sediment for the Mohr–Coulomb sliding law will be set equal to 45${}^{\circ }$. For the evaluation of the effective contact pressure, Tulaczek's et al. model (2000) will be employed.

Modelled Alpine ice-sheet extent obtained with computations using a 2 km grid (in steady state after 10 000 model years). The intensity of the blue (see colour bar) indicates the surficial ice speed in m/a. For comparison with the geomorphological reconstruction (Fig. 4) the matching rectangle is shown with a dashed black line.

[Figure omitted. See PDF]

The positive degree day (PDD) (Braithwaite, 2008) model that is implemented in PISM will be used with PDD parameters of 3 mm(${}^{\circ }$(C)${}^{-\mathrm{1}}$ d${}^{-\mathrm{1}}$) for snow and 4 mm(${}^{\circ }$(C)${}^{-\mathrm{1}}$ d${}^{-\mathrm{1}}$) for ice. The PDD model will be subjected to climate data from the WorldClim (Hijmans et al., 2005). In doing so, the same weather parameterization as in (I) is used, which separates the climate exposures of the northern and southern Alps. The climate driving is the same as that in (I): more specifically, a drop in the mean atmospheric temperature by 12 ${}^{\circ }$C is supposed, as the intensity of the precipitation reduces by 20 % on the northern side and 47 % on the southern side of the Alps. This corresponds to the parameters introduced in (I) of $\mathrm{\Delta }T=-\mathrm{12}$ ${}^{\circ }$C, $c_{\mathrm{N}}=\mathrm{0.2},c_{\mathrm{S}}=\mathrm{0.47}$. The spatial resolution of the model will be implemented with a quadratic grid of 2 km (as opposed to, 5 km in (I)). This corresponds to a 6.25-fold increase in the areal resolution. The simulations start with an ice-free topography which will be chosen to agree with today's Earth surface topography. As soon as a steady-state ice sheet geometry is reached (under time-independent climate input) the simulation is stopped and the results are analysed. The criterion for how long this computation should last is the constancy of the total ice volume (with an error of less than 3 %).

The above-mentioned parameterization will subsequently be used as a reference. Various deviations from the parameterization of this reference will be looked at with the aim of scrutinizing the effects of these alterations with regard to the build-up of the ice sheet to steady state. To discuss the role of the overestimation of the ice thickness in the Rhone Valley, our primary interest lies in the ensuing analysis of the impact of the parameter variations on the altitude of the modelled free ice surface as well as in the consequential flow dynamics of the ice. Both properties will be validated with the use of geomorphological evidence. To this end, we shall use the geographical map of Bini et al. (2009), which is again based on evidenced glacier flow boundaries. Information about the ice-flow dynamics is provided by glacial erratics, traces of striations and roches moutonnées. For the valuation of the modelled flow dynamics, the tangent to the streamlines is evaluated at all net points of the regular grid and then evaluated with the use of the locations of the origins and deposits of the found erratics.

Ice thickness

Computation in (I) indicated that the lateral extent of the Alpine ice sheet can only be adequately modelled in conjunction with a simultaneous massive build-up of the High Alps. Else, a conspicuous deviation of the ice thickness, primarily in the high Alpine valleys, is observed in comparison with the geomorphological evidence by Bini et al. (2009). This discrepancy is particularly observed in the Rhone Valley with a difference of about 800 m. In the following, several different possible causes for such deviations of the modelled ice thickness in the Upper Rhone Valley (Goms) will be looked at.

Reference simulation

Ensuing considerations are based on the previously described modelling of the Alpine LGM glaciation by use of a model climate with a temperature drop of $\mathrm{\Delta }T=-\mathrm{12}$ ${}^{\circ }$C, a reduction of the precipitation using $c_{\mathrm{N}}=\mathrm{0.2}$ in the north and $c_{\mathrm{S}}=\mathrm{0.47}$ in the southern side of the Alps and a 2 km grid. Figure 3 illustrates the corresponding spatial steady-state ice extent. In comparison with the results of Bini et al. (2009) (see Fig. 4) a conspicuous discrepancy is apparent in the thickness of the ice in the Upper Rhone Valley. Whereas the free surface of the ice is stated to be approximately 2400 m a.s.l. in the geomorphological map of Bini et al. (2009), the modelled height of the free surface of the ice is approximately 3200 m a.s.l., corresponding to a discrepancy of $\sim$ 800 m. The modelled height of the ice in this comparison is obtained from a simulation of the entire Alpine region with a quadratic grid of 5 km in length (see later Fig. 26).

Detailed image from the geomorphological reconstructed ice thickness map of the Upper Rhone Valley. Excerpt from Bini et al. (2009).

[Figure omitted. See PDF]

One recognizes increasing surficial ice velocities along the Rhone Valley. Scrutiny of the basal temperature (Fig. 5a) indicates temperate ice at the base of the Rhone Valley. The sliding velocity at the glacier bed (Fig. 5b) amounts to about 3 m a${}^{-\mathrm{1}}$ between Brig and Martigny, where its speed grows down-valley to reach more than 1 km a${}^{-\mathrm{1}}$ below Martigny. In the high mountains of the Valais and the Bernese Alps, and in the Mont Blanc Massif the ice is cold as one would expect. The modelled “equilibrium line altitude” lies in the environs of Martigny at an altitude of about 1400 m a.s.l.; this agrees quite well with the geomorphological maps of Bini et al. (2009).

(a) Distribution of the patches of temperate basal ice (blue) displayed within the relief of the glacier extent. (b) Sliding speed at the glacier bed in m/a. B is Brig, FP is Furka Pass, JJ is Jungfraujoch, M is Martigny, MB is Mont Blanc, MR is Monte Rosa

[Figure omitted. See PDF]

Influence of the model resolution

Let us study the influence of the grid resolution of the discretized model to the ice thickness. This is a valid attempt, given the approximately 800 m discrepancy obtained with Bini et al. (2009). We also note that the simulations of the entire region of the Alps with a 5 km grid in the Rhone Valley led to ice surface levels of approximately 2800 m a.s.l. To intensify the influence of the grid size, we shall reduce the 2 km resolution (Fig. 3) to the 1 km grid (Fig. 6). A comparison of these findings indicates that even this finer resolution to 1 km does not contribute to an appreciable reduction of the ice thickness. The only manifestation compared to the results with the 5 and 2 km grids is at best a somewhat more detailed pronouncement of the glacier tong that goes along with a decreased growth of the glacier length (< 5 %). It appears that no significant influence on the ice thickness is obtained by a refinement of the model resolution.

Modelled ice extent obtained for a reference run with a 1 km net in its final steady state after 10 000 model years. The blue colouring of the ice marks the intensity of the surface flow speed in m/a.

[Figure omitted. See PDF]

Influence of the precipitation and the temperature

We emphasize that the applied precipitation rate and temperature do not significantly influence the modelled extent of the height of the ice, but this holds only provided the lateral ice extent can also be obtained with the temperature and precipitation employed. This is underlined by climate studies of the entire Alpine region presented in (I) and explicitly shown in Fig. 5 in Becker et al. (2016). It shows that the climate scenarios, in which the lateral ice extent of the geomorphologic reconstruction is reached, generate free ice surfaces, which generate very similar ice heights.

Dry ice domes

The construction of the theoretical model (PISM, with the PDD climate parameterization, Braithwaite, 2008) employs a temporally constant climate with a superposed periodic annual cycle. In the PDD model a lapse rate correction of $-$6 ${}^{\circ }$C per 1000 m is used, but there is no correction of the rate of precipitation as a function of the free-surface altitude. Consequently, the precipitation remains temporally constant (corresponding to an annual constancy of the mean precipitation rate subject to monthly variation). It ignores the fact that during ice sheet build-up the free-surface height will grow, which ought to drag with it a change in the accumulation rate. Furthermore, this modelling prerequisite suppresses possible alterations of the accumulation processes that may be induced by strong winds at high altitudes. Consequently, to find out what influence this may have on the ice altitude in the Rhone Valley, the following experiment of a two-fold climate has been applied: initially the ice sheet will be built as before, but this is stopped once the lateral ice sheet extent is reached (after 12 000 model years). Subsequently, the precipitation rates above the ice domes of the Valais and Bernese Alps and in the Mont Blanc region are substantially reduced (by 20 % of the initial precipitation rates; see Fig. 7a).

(a) Distribution of the modelled differences between the steady free-surface level, obtained with the reduced precipitation rate described in the main text, and that of the reference simulation shown in blue. In the zones surrounded by the black polygons the precipitation rate is reduced to 20 % of that of the reference run; this is done once the dome regions in the Valais and Bernese Alps (VA, BA) and the Mont Blanc region (MB) are securely established. (b) Difference between the modelled ice height and that of the reference simulation.

[Figure omitted. See PDF]

Figure 7b shows the difference between the modelled ice height in steady state (after 12 000 years) and that of the reference simulation. The panel shows a decrease in the ice height from 100 to 200 m in the region between Brig and Martigny. In the remaining valley stretch, around today's region of Lake Geneva, the reduction of the ice height from that of the reference simulation is about 300 m. From now on the positions of Martigny and Brig will appear in all figure panels to easily identify the figures. The air-line distance corresponds to 40 km.

Dry Rhone Valley

A relative increase in the (lateral) advection of ice in the Rhone Valley likely contributes to the ice flow away from the domes. This could occur in view of the heavily glaciated high mountain areas at the time of the LGM. A strong build-up of the ice domes during LGM, triggered by the increase in the heights of the mountains, suggests a possible blockage of the moisture transport into the Rhone Valley. This is expressible as a reduction of the precipitation. To examine what influence a drying-out of the Rhone Valley might have for its ice height, the following simulation will serve as an estimation for the upper bound of this effect. To this end, the ice domes are built up in an initial step by the application of the reference climate (i.e. $\mathrm{\Delta }T=-\mathrm{12}$ ${}^{\circ }$C, $c_{\mathrm{N}}=\mathrm{0.2}$, $c_{\mathrm{S}}=\mathrm{0.47}$). In a second step a climate scenario is employed which reduces the precipitation rate in the region of the Rhone Valley to the 20 % level of the pre-existing precipitation rate (Fig. 8a). The reduction to the 20 % level represents an extreme case of a particularly strong desiccation. Figure 8b illustrates the deviation of the ice level once the establishment of the steady-state glaciation has been reached. The modelled ice level in the Upper Rhone Valley then lies about 200 to 300 m below the reference simulation.

Distribution of the ice level within the rectangular region of Fig. 1 for climate modelling of the glaciation of the region Valais. (a) In the indicated domain (black polygonal), i.e. domes close to this domain, the precipitation rate is reduced to 20 % of its original amount. The intensity, comprising the Rhone Valley, includes ice domains close to this domain; the precipitation rate is reduced to 20 % of its original amount. The intensity of the blue indicates the level of the ice surface in m a.s.l. (b) Deviation of the ice height from the case without reduction of the precipitation rate. The colour bar indicating this deviation in metres.

[Figure omitted. See PDF]

Simulation of the entire Würm glaciation

The preceding simulations always employed an external idealized climate, which was kept constant (but was subject to seasonal variations within an annual cycle). This led to a monotonic build-up of the ice masses until a stationary state with constant ice volume was reached. It is worth noting that in the simulations with a constant climate input, the ice volume indeed approaches a constant value but not the spatial distribution of the ice. The glacier snouts advance with unchanged frequency and they retreat, on the basis of binge–purge effects; see Mac Ayeal, 1993. It is thus reasonable to assume that the realistic climate in the Alps during the LGM was not constant but subject to natural variations, as e.g. seen in the Greenland ice core at the inception of the last ice age (North Greenland Ice Core Project members, 2004; Jouzel et al., 2007). Climate variations in the Alps during the last ice age are also manifested by a multitude of glacier advances and retreats, which have their moraine imprints in the landscape (e.g. for the Rhine glacier (see Keller and Krayss, 2005).

Deviation of the ice volume (black) and the NGRIP temperature variation (blue) plotted against time (in years before present) over the past 12 000 years.

[Figure omitted. See PDF]

It transpires that the LGM glaciation was built with several past glacier advances and retreats. Consequently, it is questionable whether the supposition of a constant climate which drives the ice sheet into a steady state is a realistic scenario. Its validity will be tested with a simulation of a realistic climate time series of the complete cycle of the last ice age. Because no such time series of the last ice age exists, a synthetic climate evolution will be generated. To this end, the variation of the signals of the oxygen isotopes of the NGRIP ice core will be applied to imitate a climate history of the Alps (North Greenland Ice Core Project members, 2004). The NGRIP oxygen isotope signal is based on time series of O${}^{\mathrm{16}}$ and O${}^{\mathrm{18}}$. These adjust to the prevailing temperature signal and favour a cooling scenario of O${}^{\mathrm{16}}$. For the ensuing considerations the variation of the signal the oxygen isotopes will be based relative to the temperature distribution $T_{\mathrm{today}}\left(x,,,y\right)$. Because the mathematical relation connecting the temperature to the concentration of the oxygen isotope ratio is non-linear, we shall employ the quadratic functional parameterization by Johnsen et al. (1995) that was also used by Seguinot et al. (2016) for the Cordilleran ice sheet evolution, $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{\Delta }{T}_{\mathrm{var}}\left(t\right)=-\mathrm{11.88}\left[\mathit{\delta },{\mathrm{O}}^{\mathrm{18}},\left(t\right),-,\mathit{\delta },{\mathrm{O}}^{\mathrm{18}},\left(\mathrm{0}\right)\right]-\mathrm{0.1925}\\ & {\displaystyle }& {\displaystyle }\left[\left(,\right(,\mathit{\delta },{\mathrm{O}}^{\mathrm{18}},\left(,t,\right),{)}^{\mathrm{2}},-,\left(,\right(,\mathit{\delta },{\mathrm{O}}^{\mathrm{18}},\left(,\mathrm{0},\right),{)}^{\mathrm{2}}\right].\end{array}$$ Moreover, we employ a lapse rate correction of $-$6 ${}^{\circ }$C per 1000 m increase in height $\mathrm{\Delta }{T}_{\mathrm{lr}}$ . Consequently, the superposition of all the temperature contributions led to the total temperature time series as given by $T\left(x,,,,y,,,,t\right)=T_{\mathrm{today}}\left(x,,,,y\right)+\mathrm{\Delta }{T}_{\mathrm{var}}\left(t\right)+\mathrm{\Delta }{T}_{\mathrm{lr}}(x,y)$ used as a climate-driving element. Figure 9 displays the generated time series of temperature reduction $\mathrm{\Delta }{T}_{\mathrm{var}}\left(t\right)$. The scaling was chosen so that the mean temperature in the time window from 32 000 to 2200 years ago corresponds to 12 ${}^{\circ }$C. In doing so, the geomorphologically reconstructed lateral extent of the ice has been adequately reproduced and the thus-modelled glaciation can be easily compared with that of the reference simulation.

Deviation of the ice height after simulation of an entire ice age, based on the oxygen isotope signal of NGRIP from the reference simulation after 28 000 years.

[Figure omitted. See PDF]

Exemplarily selected time slices of the ice build-up process (showing the free ice surface) for a simulation over the entire ice age once the maximum stage has been reached, i.e. from 30 000 years BP until today. The interval of the level lines is 100 m.

[Figure omitted. See PDF]

After the conduction of a simulation of the complete ice age from inception 120 000 years ago to the temperature optimum 28 000 years ago, the modelled ice height in the Rhone Valley is again compared with the geomorphologically reconstructed ice sheet. The chosen end of the simulation at the time of 28 000 years is based on data from NGRIP oxygen isotope signals and therefore does not agree with the instant of the LGM; in the Alps 20 000 years ago. The deviations are displayed in Fig. 10; no significant change of the ice thickness is recognized. The reason lies in a strong coupling of the ice volume with the temperature variations. The time evolution of the ice volume is compared with the temperature signal in Fig. 9. It indicates a strong volatility. This means that the glaciation must be newly built up with every new advance from near ice-free conditions, so that hardly any long ice history can be built up. Therefore, nearly no difference in the ice sheet geometry of this synthetic time-varying scenario with that of the reference simulation subjected to constant climate driving can be recognized. Examples for the responsive behaviour of the ice sheet over a time period of 120 000 years are collected in Fig. 11a–i in nine time slices of the glaciation of the Valais area. The panels show strong fluctuations of the glaciated area on both short timescales and over the entire period of 120 000 years.

Exemplarily selected modelled ice coverage time slices of the glacial disintegration process for a simulation of the entire ice age after the maximum extent, 30 000 years BP. The ice-free surface is shown with level lines at 1000 m intervals.

[Figure omitted. See PDF]

On long timescales, the glaciation of the Alps is repetitively built up and destroyed, as has already been evidenced in the time series of the ice volume in Fig. 9. Remarkably, in the ice disintegration processes the ice of large glaciers is often nearly completely melted. The computations also indicate that among similarly large Alpine glaciations the lengths of the advanced Rhone Glacier can differ. The panels b, e, h of Fig. 10 are exemplary of this fact. One recognizes that, despite the similar ice extents in the Bernese and Valais Alps, the Rhone Glacier in the lower valley is distinctly built. The reason for this can be traced back to the binge–purge effects, which after a phase of ice build-up, are associated with changes in the thermal regime and thus give rise to an enhanced sliding at the bed, which transports considerable ice mass to lower altitudes, where the ice melts. In the next cycle, the ice is then again newly built. These processes take place on short timescales and lead to strong fluctuations at the boundaries of the glaciated areas. To demonstrate the behaviour of the modelled glaciation at short timescales that is dominated by binge–purge effects, we refer to Fig. 11d, e, f. Panel d, at the time 90 000 years BP, evidences a modelled glaciation of the high-altitude regions. The Rhone Glacier has completely retreated, but grows considerably in $\sim$ 4300 years (panel e) and advances extensively into the middle of the land. Shortly thereafter (i.e. after an additional 700 years, panel f), it has again massively retreated. This climate excerpt supports the observation, expressed in connection with Fig. 9 (namely the pronounced fluctuation of the extent of the ice), which is generated by the volatility of the temperature time series. After the maximum extent of the modelled ice sheet that occurred about 30 000 years BP, the final glacier retreat into the Holocene commences (Fig. 12).

Difference in the modelled ice height for the simulations after the variation of the seasonal air temperature parameter $\widehat{T}=\mathit{\{}-\mathrm{5}{}^{\circ }\mathrm{C}$ (a); $-\mathrm{2}$ ${}^{\circ }$C (b); +2 ${}^{\circ }$C (c); +5 ${}^{\circ }$C (d)}. The lateral extent of the ice is considerably reduced for $\widehat{T}=+\mathrm{5}$ ${}^{\circ }$C (blue line, panel d) compared to the geomorphological reconstruction of Schlüchter et al. (2009), in which the entire area of panel (d) is covered by ice.

[Figure omitted. See PDF]

Effects of the seasonal dependence of the climate

An additional cause for the high ice altitudes in the Rhone Valley could also be a larger seasonal dependence of the air temperature during the last ice age. This thesis is supported by Coleopteran findings from the time of the marine isotope stage 3 about 57 000 years ago. These findings are likely due to high summer temperatures in July as well as extraordinarily low temperatures of $-$15 to $-$7 ${}^{\circ }$C in January–February. If one compares these values with air temperatures of today (Becker, 2017, Fig. 1.2), these data correspond to an increase in the amplitude of the seasonal temperature variation of approximately 5 ${}^{\circ }$C. Even though no explicit evidence is known for the existence of such a mean air temperature increase for the time at the LGM, we shall nevertheless consider the hypothesis here. To this end, the seasonal temperature amplitude will be varied. This amplitude will be increased and reduced respectively to fathom the corresponding impacts to the ice sheet extent and depth distribution. Consequently, one obtains the following dependence of the air temperature on time $t$ (in months) and on spatial coordinates, $x,y$, $$\begin{array}{cc}{\displaystyle }& {\displaystyle }T\left(x,,,,y,,,,t\right)=T_{\mathrm{today}}\left(x,,,,y\right)+\mathrm{\Delta }T+\mathrm{\Delta }{T}_{\mathrm{lr}}\left(x,,,,y\right)\\ & {\displaystyle }& {\displaystyle }+\mathrm{\Delta }{T}_{\mathrm{season}}\left(t\right),\end{array}$$ in which $\mathrm{\Delta }T$ is the temperature correction (by $-$12 ${}^{\circ }$C, relative to the reference simulation, discussed in Sect. 4.1). Furthermore, a linear cooling with altitude is applied with a lapse rate of $-$6 ${}^{\circ }$C per 1000 m, expressed above by $\mathrm{\Delta }{T}_{\mathrm{lr}}$. In addition, a seasonal correction of the temperature, denoted by $\mathrm{\Delta }{T}_{\mathrm{season}}\left(t\right)$, is introduced. This correction is given by $$\mathrm{\Delta }{T}_{\mathrm{season}}\left(t\right)=\widehat{T}\sin \left(t,-,\mathrm{4}\right)\left(\text{time}t\text{in months}\right),$$ in which $\widehat{T}$ is the amplitude of this seasonal variation. It will be chosen as $\widehat{T}=\left\{-,\mathrm{5},,,,-,\mathrm{2},,,,+,\mathrm{2},,,,+,\mathrm{5}\right\}{}^{\circ }$C. To obtain the maximum amplitude in July, the phase of the sine function was shifted by 4 months. Figure 13 displays the difference in the modelled ice level when the dependence of the air temperature deviates from that of the reference simulation by $\widehat{T}=\left\{-,\mathrm{5},,,,-,\mathrm{2},,,,+,\mathrm{2},,,,+,\mathrm{5}\right\}{}^{\circ }$C. Thereby the lateral ice extent decreases conspicuously with an increase in the seasonal dependence of the air temperature (Fig. 13c, d). A considerably larger deviation is evidenced at low-lying land, where the ice height is strongly reduced in comparison to the surface level of the reference simulation. Along with this, the ice thickness in the Rhone Valley remains essentially unchanged. Only once the amplitude $\widehat{T}$ is increased by $+$5 ${}^{\circ }$C (Fig. 13d) a small decrease of 100 to 200 m of the glacial height is observed in the region around Martigny; this is paired with a considerable reduction of the horizontal extent of the ice.

(a) Glacier bed topography after removal of the sediments and/or basal rocks. (b) Deviation of the ice height from that of the reference simulation after removal of the sediments in the Rhone Valley.

[Figure omitted. See PDF]

The counter-evidence is shown by the results obtained with a reduced seasonal dependence of $\widehat{T}$ by $-$5 and $-$2 ${}^{\circ }$C (Fig. 13a, b). In these cases, an increase in the ice coverage is manifested in the Alpine foreland north of Martigny. The ice thickness in the Rhone Valley remains essentially unchanged. So, we may conclude from these findings that the lateral ice extent is considerably reduced when the seasonal variation of the air temperature increases, whereas it strongly increases when the seasonal variation of the air temperature decreases. Thus, an increase in the seasonal temperature variation worsens the results, because the influence of the lateral ice extent is of greater impact than the reduction of the ice height.

Sedimentation

Apart from climatic causes, we also wish to scrutinize topographic causes of the studied discrepancies of the ice sheet heights in the Rhone Valley. For instance, the implemented present glacier bed topography could considerably deviate from the bed topography shortly before the LGM. A possible cause could be the deposition of sediments at the valley base. Dürst Stucki and Schlunegger (2013) constructed a topography of the rocks of the central Alps in a stratigraphic exploitation of more than 40 000 boreholes. In their report, the sediment filling of the Rhone Valley was estimated to be larger than 1000 m thick. This means that the deepest point of the Rhone Valley without the sediment filling must have been approximately 400 m below today's sea level. How much the valley was exactly filled with sediments at the time before the LGM is not known. This is why today's DEM was used as a glacier and basal rock bed at the time of the LGM. To estimate the influence of the sedimentation on the ice height, we shall now (as an extreme case) artificially remove the complete sediment bed in the Rhone Valley. Figure 14a shows the employed sediment-cleared topography of the glacier bed.

The modelled steady ice height so obtained agrees in all details with the reference simulation. Figure 14b indicates that the removal of the sediments does not exert an appreciable influence of the ice height in steady state. This means that the ice thickness is increased by the amount of the removed sediments; yet, the level of the free surface remains surprisingly robust against a removal of sediments.

Grid-point densification at the knee of the Rhone River at Martigny

The best grid resolution with which the simulations can be performed in reasonable CPU time is about 1 km. This choice implies that at the abrupt 90${}^{\circ }$ bend of the Rhone knee at today's position of Martigny only about 10 grid points are at our disposal for modelling the transport of the ice out of the Rhone Valley in the direction of Lake Geneva. This suggests that the excessively small number of grid points causes artificial clogging of the ice transport. To study the influence of the small number of grid points in the region of Martigny, we apply a glacial bed topography in a first simulation, in which the western flank of the Dents de Morcles was artificially removed. So, a considerable widening of the Rhone Valley close to Martigny is induced. The through-flow cross-section area is roughly doubled, which furthers the ice transport through the valley toward Lake Geneva. In a second simulation, a glacier bed topography is applied, and in addition to the Dents de Morcles, the whole of the Muveran mountains are artificially removed. The topographies of both simulations are indicated in Fig. 16. The second simulation represents an upper extreme case in which the large widening of the Rhone Valley enlarges the consequential increase in the through-flow cross section by a factor of 5.

Figure 15a, b display at steady state the deviations of the ice height from the reference simulation. Panel a shows that the artificial removal of the western flank of the Dents de Morcles reduces the steady ice height in the Upper Rhone Valley relative to the height of the reference simulation by about 30 m. For the case of the removal of the complete Muveran Mountain, the modelled ice height relative to that of the reference simulation decreases by about 100 to 200 m.

Deviation of the modelled ice height from the reference simulation with the removal of the western flank of the Dents de Morcles (a) and of the whole of the Muveran mountains (b).

[Figure omitted. See PDF]

Representation of the artificially removed rocky regions of the Dents de Morcles (light blue) and the Muveran mountains (dark blue) ©swisstopo.

[Figure omitted. See PDF]

Deformability of the ice

In Dye 3 ice-core analyses from the GRIP ice core, Shoji and Langway (1987) found an increased deformability of the ice that is larger by a factor of approximately 3.5 in a layer that is attributed to the last glaciation (Würm). Thus, the fluidity of Glen's power law is increased by an enhancement factor $E=\mathrm{3.0}$. The increased deformability is caused by impurities in the ice due to an enlarged concentration of dust (Hutter and Vulliet, 1985). In general, the increased concentration of dust in the Greenland ice is attributed to enhanced global wind speeds during the LGM, which transported the dust over large distances. For ice of the last glaciation the enhancement factor is on average about 2.5 times larger than for Holocene ice (Coffey and Paterson, 2010). To test which influence the deformability of the ice exerts on the ice thickness in the Upper Rhone Valley, the modelling is extended by the introduction of an enhancement factor into Glen's flow law. The panels in Fig. 17 display the deviation ice heights (for steady-state conditions) and maintained reference simulations when the enhancement factors $E=\mathit{\{}\mathrm{3},\mathrm{10},\mathrm{20}\mathit{\}}$ are used. Computations indicate $E=\mathrm{3}$, measured in a Greenland ice core, and a reduction of the modelled ice height in the Upper Rhone Valley by approximately 100 m relative to the modelled height for the standard simulation with $E=\mathrm{1}$. Use of an enhancement factor of $E=\mathrm{10}$ reduces the ice height by 300 to 400 m. Only when $E=\mathrm{20}$ is used, the height of the modelled ice sheet ends with a height reduction of 500 m compared to the standard simulation of the geomorphological reconstructed ice height. However, this value of $E$ is unrealistically large.

Validation using geomorphological evidence

To evaluate the level of the free ice surface from geomorphological evidence, we used the differences in the digitized DEM of the modelled ice height and the DEM of today's topography of the Alps. With this, the nunataks that surmount the free ice surface can be identified. These were compared with the nunataks along the boundaries of the striations listed by Bini et al. (2009). For better comparability, the modelled nunataks were directly projected onto the geomorphologic map. Results are shown in Fig. 18.

Deviation of the ice height for the reference simulation using three different enhancement factors and the standard reference simulation ($E=\mathrm{1}$): $E=\mathrm{3}$ (a), $E=\mathrm{10}$ (b), $E=\mathrm{20}$ (c).

[Figure omitted. See PDF]

Basal area of the modelled nunataks (red) superimposed on the geomorphologically reconstructed map of the ice extent during the LGM due to Schlüchter et al. (2009). The black dashed boundaries of the rectangular regions are represented in detail in Figs. 19 and 20. The upper (larger) region shows the Bernese Alps in Fig. 20; the lower (smaller) region shows the Valais Alps in Fig. 19.

[Figure omitted. See PDF]

Enlarged and more detailed map of the Valais Alps, close-up of Fig. 18.

[Figure omitted. See PDF]

Enlarged and more detailed map of the Bernese Alps, close-up of Fig. 18.

[Figure omitted. See PDF]

If one ignores the adulterated regions along the boundaries, the figure shows a smaller basal area of the modelled nunataks compared to those of the geomorphological map in the background, caused by the higher modelled ice surface. Owing to this, some of the nunataks of the map cannot be reproduced in the model (e.g. Bella Tolla). If one inspects and compares the details of the Figs. 19 and 20 of the Valais and Bernese Alps, one recognizes that, particularly in the Valais Alps, only a few mountain peaks and mountain chains surmount the ice. A better agreement prevails in the Bernese Alps, where the modelled ice heights in the northern part of the modelled region agree better with the corresponding geomorphologically reconstructed ice heights than in the southern part.

Discussion

It is concluded that the modelled glaciation in the Alps (Fig. 3) during the LGM hints at an ice-dome-like topography of the glaciation – also in the western Alps – rather than at an ice-flow-net-like topography. In a direct comparison with the geomorphologically reconstructed glacial extent from the land map of Bini et al. (2009), it is noticeable that the Rhone ice dome is less pronounced and may be considered to be a kind of a foreland of the Bernese ice dome. Similarly, the Mattertal ice field of the land map is also a less independent ice mass in the modelled topography and is instead covered by an ice dome with its centre above the Monte Rosa region. Both ice domes strongly dominate the topography of the glaciation of the Upper Rhone Valley, which is accompanied by a distinct ice divide between the Bernese and Valais ice domes.

Influence of the various changes in the model ice heights in the Upper Rhone Valley.

Model change	Reduction of the ice height in the Upper Rhone Valley
Refinement of grid resolution by 5, 2, 1 km in length	–
Ice domes subject to dried climate	100–200 m
Dry land climate in the Rhone Valley	200–300 m
Time series of the driving climate of the last ice age	–
Increase in the seasonal temperature variation	–
Sedimentation	–
Removal of the throat of the Martigny knee	100–200 m
Flow enhancement by dust in the ice ($E=\mathrm{3}$ versus $E=\mathrm{1}$)	100 m

As a consequence of this ice-dome-like glaciation, a deviation of the ice height and ice thickness between the computed and geomorphological reconstructions can be recognized. The detailed study reported here focuses on the modelled ice height. In the Upper Rhone Valley, it could be quantified as roughly 800 m above the geomorphologically reconstructed ice height of the land map of Bini et al. (2009). Various alternative assumptions regarding climate, topography, grid size of the discretized model and rheology of the ice were implemented with regard to the Upper Rhone Valley to explain the discrepancy of the ice height predicted by the PISM-based computations and the inferences by Bini et al. (2009) and to reduce them. The simulations, however, indicated no significant differences of the modelled ice height, even when results were produced with the grid sizes of 5, 2 and 1 km. (Sect. 4.2).

That the increased grid resolution alone was not able to explain a height difference of 800 m between our model results and the geomorphological inferences was enforced by computations described in Sect. 4.9. In these simulations, the mountains at the knee of the Rhone River at the height of today's Martigny were artificially removed in order to enlarge the number of grid points in the throat of the valley and to lower its blocking power. The artificial removal of the Dents de Morcles corresponds roughly to a doubling of the ice-flow cross section at the knee of the Rhone Valley; however, this widening generated no more than a 30 m ice height reduction. By removing the complete Muveran mountains, it was possible to reduce the ice height by 100 to 200 m. Consequently, it was shown that the relatively small grid point density in the environs of this narrow pass did not provide a sufficient explanation for the 800 m ice height difference.

Furthermore, it became evident that a spatially constant adjustment of the precipitation rate and the temperature variation did not result in significant changes in the ice height in the Upper Rhone Valley. Rather, a reduction of the precipitation rate above the ice domes evidenced implications for the ice heights in the valley (Sect. 4.4). For example, a desiccation of the Mont Blanc Massif and of the Bernese or Valais Alps led to a reduction of the ice height by 100 to 200 m in the Upper Rhone Valley. Analogously, the implications to the ice height of a possibly emerging desiccated inland climate in the Rhone Valley were studied in Sect. 4.5. With these, it could be shown that the desiccation of the Rhone Valley after the LGM extension led to a reduction of the ice height by 300 m.

All above discussed simulations were based on an annual cycle with a temporally constant climate; this was reflected in a monotonous build-up of the ice masses. In order to study the influence of a constant climate on the build-up process and the overestimation of the ice thickness, a synthetic LGM climate was employed that was imprinted with additional natural seasonal variations. In Sect. 4.6 a complete ice age was modelled that was based on a time series of variations of the NGRIP oxygen isotope O${}^{\left(\mathrm{16}\right)}$ $/$ O${}^{\left(\mathrm{18}\right)}$-ratio. The employed temperature time series was scaled so that the lateral ice extent at the time of the climate optimum corresponded to that of the geomorphologically reconstructed ice height distribution and could be compared with the foregoing modelling attempts. The results of the simulations made it clear that the inferred ice heights in the Upper Rhone Valley are the same as those of the reference simulation (of constant climate). Consequently, the discrepancy with the geomorphologically reconstructed ice height also remained unchanged when a complete variable climate time series of the past 120 000 years was used.

As a further climatic cause for an ice height configuration change, a conspicuous seasonal dependence of the air temperature in the Rhone Valley was considered. Yet, it was shown in Sect. 4.7 that the seasonal variation of the driving climate affects the ice height in the Rhone Valley only minimally. Instead, it shows a large influence of the lateral ice extent. This can be attributed to the high summer temperatures, which induce ice melting at low altitudes of the Alpine foreland.

Moreover, the results obtained in Sect. 4.8 have shown that the modelled ice height is equally robust against changes in the basal sediment filling in the Rhone Valley. It was evidenced that a reduction of the basal sediment filling leads to an ice thickness increase comparable to the removed sediment mass, which leaves the ice surface level practically unchanged. By contrast, a change of the deformability of the ice by an alteration of the enhancement factor is paralleled by a change of the ice thickness. If the enhancement factor is set equal to $E=\mathrm{3}$, as supported by specimens of the Dye 3 ice core (see Sect. 4.10), the modelled ice height drops by roughly 100 m when compared with the reference simulation. Unfortunately, for the free-surface ice level in the Rhone Valley none of the discussed hypotheses or influence factors suffice by themselves to explain the extent of this drop. However, Table 1 suggests that at least a combination of local climate changes during the LGM and beyond, as well as an inadequate grid resolution of the model in the vicinity of Martigny combined with an enhanced deformability of the LGM ice height predictions, may, perhaps, partly explain the discrepancy of the ice height predictions of the model versus the geomorphologically reconstructed alternative.

Ice dynamics

In the preceding section, our focus was on how the deviations between the modelled ice height and its geomorphological reconstruction could be reduced. In order to compare the model results with the aid of further geomorphological evidence and to scrutinize the consequences of the ice height onto the local flow regime of the ice age glaciers, here we shall consider the flow field of the glacial ice sheet in the Valais region in detail. We shall construct flow lines for the reference simulations so far employed in order to visualize the flow field. Moreover, we shall consider the flow field during the ice build-up process. In addition, the influence of the employed distribution of the precipitation on the flow field will also be studied. In particular, the region of the Simplon Pass and the Monte Leone–Blinnenhorn chain (as displayed in Fig. 2) will be looked at. This will be open to questions regarding the work of Bini et al. (2009) and Kelly et al. (2004). These authors formulated the north–south transfluence of ice during the LGM. By this we mean that a discharge of ice across an ice-sheet takes place, even though under ice-free conditions such a transfluence is not possible for water. The mathematical prerequisites of the so-called shallow ice approximation are such that ice flow in an ice sheet is always down the steepest descent of the free-surface topography. So, if an ice divide crosses a north–south-oriented water shed boundary, ice will flow from it to the east and the west, in contrast to the functioning of a water shed. The ice flow will be evaluated with the aid of motion of the erratic boulders.

Free-surface ice-flow lines (white) overlain on the topography of the glacier bed after a simulated ice build-up into a quasi-steady state (after 12 000 years). The thick white boundary line shows the extent of the glacier that was reached.

[Figure omitted. See PDF]

Nine time slices of the ice build-up evolution (free-surface level) over 12 000 years. The contours underlain in blue indicate the free-surface height in steps of 500 m intervals. The white colouring corresponds to a height of more than 4000 m. The stream (flow) lines drawn as red solid lines show the flow directions of the ice (all over the ice depth). The dashed thick red lines mark the resulting ice divide.

[Figure omitted. See PDF]

Time slices of the ice build-up evolution of the ice sheet under present driving precipitation scenarios. The contours shown in blue mimic the free-surface topography in steps of 500 m. Thereby, the white shading corresponds to a free-surface height of more than 4000 m. The red solid lines show the stream (flow) lines. The dashed, thick red line marks the location of the ice divide. The times in the inset boxes denote the time slices after inception of the ice build-up.

[Figure omitted. See PDF]

Ice-flow dynamics under steady-state conditions

Figure 21 displays computed streamlines (see Methodology); they are underlain by the topography of the glacier bed. The ice of the Rhone Valley is discharged into two regions: whereas the major part of the ice flow is essentially downstream from Brig, following the Rhone Valley, the glacial ice above Brig flows toward the Goms region and further south over the Monte Leone–Blinnenhorn chain. This means that part of the ice, contrary to the present discharge of the Rhone River, does not follow the topography of the Rhone Valley, but follows the ice transfluence over the mountain chain lying orographically left, against which it seems to be “pressed”. More explicitly, ice from the Goms area is transported into the Italian Val d'Ossola. In terms of the SIA ice approximation, this means that the ice dome or divide at the LGM lies to the south of the topographic north–south watershed.

The separation into two ice discharges is thus marked by an ice divide that delimits the boundary of the two regions of the ice discharge and does not agree with the watershed boundary. Its location connects the two ice domes of the Bernese and Valais Alps (Fig. 3). The Bernese and Valais ice streams discharging into the Rhone Valley form a diffluence here, which separates the Rhone Valley into two parts. The ice in the lower part of the Rhone Valley flows in a downstream direction, whereas the ice on the other side of the ice divide flows over the Monte Leone–Blinnenhorn chain, thus leaving the desiccation system of the Rhone River. Note that this is a strict inference of the SIA employed in PISM and can only be corroborated by data, e.g. by erratics. If their positions do not provide support for this split into two discharge regions, the results of the PISM computations will have to be replaced by software integrating of the full Stokes equations.

Development of the ice surface

In the last subsection, the flow field of the Alpine glaciation in the Upper Rhone Valley was looked at under steady-state conditions. In order to understand how the flow regime of the transfluence from the Goms into the Val d'Ossola can be formed, we shall now study the temporal development of the surficial flow lines during the ice-sheet build-up process. Figure 22 shows the development of the ice build-up, starting with initial ice-free conditions for which the flow over the glacier bed topography is prescribed up to the established transfluence toward the end of the ice build-up. The individual time slices demonstrate how the concurrent spatial displacements of the ice divide develop. In the simulated run the same reference parameterization of the climate will be used (see Sect. 4.1: temperature drop of 12 ${}^{\circ }$C, reduction of the precipitation in the north reduced to a level of 20 % and in the south to 47 %). The simulation will be stopped as before as soon as a stationary state of the glaciation has been reached.

The panels in Fig. 22 show that during the time span between the inception of the simulation and 1000 model years afterwards, the ice divide follows a trace from the Monte Rosa Massif to the Furka Pass, analogously to the water shed today. Panel a shows an ice dome near the Furka Pass, which corresponds to the ice dome in the land map of Bini et al. (2009). Already at this time the ice height corresponds to that of the Rhone ice dome on the land map of Bini et al. (2009), approximately at 2000 m. At this time slice, no transfluence across the Simplon Pass has yet been established, despite the fact that the geomorphologically reconstructed ice thickness had already been reached. After approximately 1300 model years (panel b), the shape of the Rhone ice dome mutates to a plateau glacier in the Upper Goms, southwest of the Furka Pass. Then, several ice domes are formed in the Bernese and Valais Alps. This plateau–glacier grows further with ongoing glaciation and shifts its centre after 1500 model years to the central Upper Rhone Valley (panel c). Consequently, the ice divide rotates further anticlockwise in the direction of the north–south axis when a connection between the Goms-Plateau glaciers is formed concurrently with the growth of the ice dome in the Bernese Alps. After about 1700 model years (panel d) the plateau glacier has merged with the Bernese ice dome. Therewith, for the first time the ice divide lies on a north–south-oriented axis between the Bernese and Valais ice domes. The streamlines now cross the Monte Leone–Blinnenhorn chain and the region of the Simplon Pass; this marks the inception of the ice transport from the Upper Rhone Valley into the water shed system of the Toce river. In the following years, the axis of the ice divide is even more north–south-oriented and finally stabilizes itself in this orientation. The ice divide is now approximately positioned in the region of the town Brig (panel e) and so separates the discharge regime of the Upper and Lower Rhone Valley. More specifically, the ice discharges into the lower, western, part of the Rhone Valley toward Lake Geneva, whereas in the upper eastern part the ice moves southward and discharges into the Val d'Ossola. Therewith, a switch of the location of the ice divide is induced from a topography-following flow configuration to a relief-independent ice divide formation. The final orientation of the Valais ice divide is almost orthogonal to the initial ice divide that appears to be topographically constrained.

Sensitivity of the ice surface topography to the climate

In this subsection, we study, how strongly the distribution of the precipitation influences the ice surface level and the flow behaviour of the ice. To this end, a further simulation is performed; it essentially reproduces the scenario of the last simulation but is not modified; i.e. the north–south separation of the precipitation pattern of the reference model of Sect. 4.1 is not used, rather the precipitation pattern of today is used. A reduction of the temperature of 12 ${}^{\circ }$C (relative to today's temperature) as well as a reduction of the precipitation rate by 40 % (of today's precipitation rate) in the northern and southern sides of the Alps will be employed in contrast to the parameterization of the climate before. With this precipitation pattern the LGM ice extent could not be reproduced. Here, we shall now study the impact that extreme precipitation has on the emerging ice sheet.

Locations of the origin of erratics (small green circles) found in the Bernese Lake District (large green circles), as well as possible trajectories of their transport routes (green dashed lines), which follow the computed flow lines (red) of the modelled ice dynamics. One recognizes the transport into the Rhone Valley. How the erratics, however, reach from their ends of the trajectories close to today's Lake Geneva to their positions in the Bernese Lake District, is not part of the present paper. The geomorphologically reconstructed flow lines in the region of the Grimsel (after Wirsig et al., 2016, yellow dashed line). The ice surface level is indicated by the blue background (shading in 500 m intervals; white shading means a height above 4000 m a.s.l.).

[Figure omitted. See PDF]

In Fig. 23 the resulting temporal development of the flow lines are displayed. Generally, a very similar ice build-up process to the previous case is generated, in which two north–south-separated climate zones were employed. One recognizes that conspicuous steady transfluences appear in the Monte Leone–Blinnenhorn chain region, as well as in the neighbouring crossings. Moreover, transfluences in the region of the Simplon do not occur, not even in the time slice at 3000 years (quite contrary to the cols of the Monte Leone–Blinnenhorn chain). At this time slice, the ice height upstream of Brig already overtops that of the geomorphological reconstruction (compare with Fig. 4, Bini et al., 2009). Even after 10 000 years, when the ice sheet shape has long stabilized, no definite transfluence across the Simplon Pass can be identified. Rather, the Valais ice divide follows a distinct curve from the Monte Rosa Massif, exactly across the peak of the Simplon Pass to the Bernese Alps. This also represents a significant difference from the results obtained with the reference modelling. Contrary to the latter, in which steady conditions are reached slightly after 3000 years, steady conditions are reached in this case only after 7000 years.

Nevertheless, when clearly different climate drivers are applied, transfluences from the Upper Rhone Valley into the system of the Toce can still develop. These other scenarios also establish an ice divide across the Rhone Valley. The ice divide is in these cases located within the region of the Simplon Pass, so that a possible transfluence cannot be clearly identified.

Validation of the modelled flow lines

To validate the modelled flow field with geomorphological evidence, particular erratic boulders can be used, the original locations of which could uniquely be determined. Figure 24 shows streamlines in the entire modelled region for stationary conditions of the reference simulation. When computing these, one had to pay attention to the fact that the deviations of the computed flow lines from the correct ones grew toward the boundaries of the computed domain. The figure also shows the points of erratics, as listed in Table 2.

List of erratics that were used for the validation of the computationally determined flow lines.

Location of origin	Rock description	Location found
Mont Blanc Massif	Mont Blanc granite	Bözingenberg near Biel (BE), Walliswil (BE)
Vallorcine	Vallorcine conglomerate	Finsterhennen (BE), Lyss (BE)
Findelen Glacier	Serpentine	Hindelbank (BE)
Allalin, Saas-Fee	Allalin Gabbro	Finsterhennen (BE)

The locations of the origin of all these erratics are in the area of the larger Mont Blanc and Monte Rosa regions. The erratics were all found in the Bernese Lake District, but the origins of the Mont Blanc granite are all on the northern part of the Mont Blanc Massif; they were likely transported northward from their origins. They all lie in the drainage area of the Rhone and were transported toward Lake Geneva. This is in good agreement with the modelled flow lines.

Even if no agreement with modelled flow lines and erratic trajectories could be found in the Lake Geneva region and in its continuation (one reason for this may be boundary effects), it is nevertheless evidenced that all erratic boulders lie in the catchment area of the Rhone. This is a prerequisite for the erratics to be deposited in the Solothurn Lobe of the Bernese Lake District (compare this with the reconstruction of the Solothurn Lobe in Bini et al., 2009). We also mention that the flow lines in the region of Oberhasli show good agreement with the reconstruction of the discharge there. This transport follows the course of the Aare river from the Grimsel Pass to the present location of Meiringen (Wirsig et al., 2016).

Detail of Fig. 2 of Kelly et al. (2004). Positions of geomorphological evidence indicating the vertical LGM ice extent is shown. For the case that a piece of evidence discloses information about the flow direction of the ice, it is drawn with a corresponding arrow, else with a point. The region discussed in the text is the rectangle displayed in red of the Monte Leone–Blinnenhorn chain. The figure shows that only very little evidence was found (red arrows); these, however, hint at a southward flow. The LGM evidence at the Simplon Pass is also shown.

[Figure omitted. See PDF]

Modelled steady-state extent of the Alpine ice sheet during LGM for the indicated temperature drops $\mathrm{\Delta }T$ and precipitation correction factors, $c_{\mathrm{N}}$ and $c_{\mathrm{S}}$, for the northern and southern Alps. The geomorphologically reconstructed extent of the ice is indicated with the red line. The level lines correspond to 1000, 2000 and 5000 m; from Becker et al. (2016).

[Figure omitted. See PDF]

Discussion

It was shown in this study that the free-surface topography of the modelled last glaciation in the Upper Rhone Valley is connected to the flow dynamics, taking place in the existing topography and the climatological external input. The model results hint at the presence of an ice divide between the Valais and Bernese Alps, which generates north–south transfluences from the Rhone Valley into the valley of the river Toce. These results show that (i) ice in the region of the Simplon Pass was transported from the catchment system of the Rhone to that of the Toce, and (ii) transfluences are being formed in the neighbouring transitions of the Monte Leone–Blinnenhorn Chain (e.g. Albrun Pass, Chrigalp Pass). In these cases too, ice from the Rhone Valley is transmitted to the Val d'Ossola. Whereas the transfluence at the time of the LGM has already been discussed in geomorphological articles, e.g. Bini et al. (2009) and Kelly et al. (2004), the transfluences in the transitions of the Monte Leone–Blinnenhorn chain have hardly been discussed in earlier articles. In Arn (1998) roches moutonnées are described, which were found at the Albrun Pass; here too a transfluence of ice from the Rhone into the Toce catchment area could be identified (Kelly et al., 2004). Moreover, evidence is described at the Geispfad, a transgression close to the Albrun Pass (Fig. 25). Whereas in this article only the flow directions are given, Arn (1998) points out that the mountain slopes of the Albrun Pass (2409 m a.s.l.) above 2600–2700 m a.s.l are more conspicuously weathered than elsewhere. He infers from this an ice depth of 200 to 300 m at the Albrun Pass. In addition, in the present simulations, a robust glaciation of the above passes is obtained. The computational results therefore agree very well with the mentioned geomorphological evidence. The use of the deduced LGM climate in previous subsections, which is characterized by strong precipitation in the southern Alps, leads to the formation of the Valais ice divide between the Bernese and Valais Alps at the position of Brig. In general, the formation of an ice divide between the Valais and Bernese Alps is surprisingly robust against the given climate. This was scrutinized by a study of the stability of the results relative to the precipitation scheme that is dominated by southern flows. This was done in a separate simulation, in which today's precipitation scheme, which is dominated by flows from the north, was employed. An ice divide is also formed for this climate; it gives rise to related transfluences from the Rhone Valley into the catchment area of the Toce river. However, in this case the Valais ice divide moves so far eastward that no transfluence at the Simplon Pass is formed. This contradicts the existence of striations found at the Simplon Pass by Kelly et al. (2004), which definitely hint at a transfluence. Therefore, on the one hand the robustness of the ice divide is demonstrated and, on the other hand, the transfluence is corroborated. This is also indirectly confirmed by the modelled LGM precipitation scheme, which is strongly imprinted by flows from the south, since the latter is also necessary to reconstruct the Simplon transfluence.

In summary, the position the Valais ice divide depends on the dominance of the precipitation on the southern side of the Alps. The more these prevail, the more the ice divide will be shifted to the west. The orientation of this divide is robust against the employed precipitation scheme; it is orthogonally oriented along the underlying Rhone Valley. This emphasizes and explains the alteration of the discharge regime driving the ice build-up process. During this process, an alteration of the ice build-up flow, governed by the mountain chain (corresponding to today's watershed), takes place toward an ice-cap-dominated flow over and across the underlying mountainous topography of the underlying mountain structure.

Moreover, the modelled flow of the Alpine glaciation was compared with locations of the origin and deposition of erratics. A satisfactory agreement was thereby demonstrated between model-induced results and geomorphological evidence. The locations of deposition of the considered erratics lie in the region of the Bernese Lake District, a fact corroborating the modelled ice flow northward. Yet, no information is known about erratics south of the passes, which could be attributed to the Upper Rhone Valley, to validate the modelled southward flow across the Simplon Pass or the passes of the Monte Leone–Blinnenhorn chain. A different reason could be that PISM, with the incorporated SIA would not be applicable for this case, as was already hinted at on several occasions.

For the height of the ice mass, discussed in Sect. 4.12, the exploitation of the ice-flow field also delivers a contribution. It supports the fact that, for an ice height in Goms of only 2000 m, as it is represented in the geomorphologically reconstructed land map, no ice will flow over the Simplon Pass from the Rhone into the Toce catchment region; it will rather follow the Rhone Valley. Only conspicuously higher free surfaces over the Simplon Pass are likely to occur. Therefore an ice-dynamical discrepancy with the geomorphological land map of Bini et al. (2009) has been discovered.

Summary and conclusions

In this work, detailed simulations of the glaciation during LGM in today's Upper Valais and the neighbouring regions have been performed. A numerical ice sheet model called PISM was employed, which is based on the SIA with incorporated Mohr–Coulomb-type basal flow properties. This model was subject to a LGM precipitation model, which is characterized by enhanced precipitation rates on the southern side of the Alps. This model shows good reproduction of the Alpine ice sheet. Unfortunately, a conspicuous discrepancy of the modelled level of the free ice surface with the geomorphological reconstruction, observed in the considered Upper Rhone Valley, could not satisfactorily be reduced. For this reason, various alterations of the model were implemented with the aim of reducing this discrepancy of $\sim$ 800 m. The discrepancy could neither be explained nor adequately reduced, because it could not uniquely be associated with a mathematical or physical cause. A combination of several alterations of the model, e.g. a refinement of the grid size, the incorporation of a snow model and softer ice and a better climate model, might well reduce the discrepancy. A possible way to improve the mathematical–physical model could be replacing PISM with the Stokes approximation; yet CPU times with these are far too large with today's computer efficiencies. Moreover, scrutinizing and questioning the inferences obtained by the geomorphological approach is equally advisable, as it is simply another approach using indirect evidence.

The model describing the flow dynamics of the ice masses in this Swiss region has shown that (subject to the admissibility of the SIA employed in PISM) the model tends to the formation of an ice divide between the ice domes of the Bernese and Valais Alps at and in the surroundings of today's Brig. This ice divide is at high glaciation roughly perpendicular to the underlain Rhone Valley, which results in two different ice discharge regimes. Whereas the ice, downstream of Brig, follows grossly along and downstream of the Rhone Valley to today's Lake Geneva, the ice upstream of Brig flows over and across the Simplon Pass, as well as across the neighbouring passes of the Monte Leone–Blinnenhorn chain into the catchment area of the Toce river. The results presented herein can thus not only corroborate the geomorphological studies, which hint at an influence from the north and south, in the region of the Simplon Pass and Albrun Pass, but also at the transfluence with the presence of a Valais ice divide between the Bernese and Valais Alps.

The raw data are tons of GB and not readable without a detailed introduction of the authors. But the prepared data are accessible in the figures of this paper as well as in more detail in the PhD thesis (Becker, 2017).

The authors declare that they have no conflict of interest.

Acknowledgements

The research profited from the financial support of D-BAUG, ETH Zürich. Support was also provided by Julien Seguinot and members of the Glaciology Section at VAW. Thanks are also expressed to the PISM authors for help with their software. Edited by: Martin Hoelzle Reviewed by: two anonymous referees