Last Glacial Maximum precipitation pattern in the Alps inferred from glacier modelling

During the Last Glacial Maximum (LGM), glaciers in the Alps reached a maximum extent, and broad sections of the foreland were covered by ice. In this study, we simulated the alpine ice cap using a glacier flow model to constrain the prevailing precipitation pattern with a geomorphological reconstruction of ice extent. For this purpose we forced the model using different temperature cooling and precipitation reduction factors. The use of the present-day precipitation pattern led to a systematic overestimation of the ice cover on the northern part of the Alps relative to the southern part. To reproduce the LGM ice cap, a more severe decrease in precipitation in the north than in the south was required. This result supports a southwesterly advection of atmospheric moisture to the Alps, sustained by a southward shift of the North Atlantic storm track during the LGM.


Introduction
In the mid-nineteenth century, systematic studies of glacial features (e.g.Agassiz, 1838;de Charpentier, 1841) yielded the formulation of the ice age theory.Although it was long unclear whether there had been a single or multiple glaciations, this controversy ended with the systematic classification of the glaciofluvial terraces of the northern forelands, indicating that there had been at least four glaciations in the Alps.However, the landscape imprint of previous glaciations has been largely overprinted by the more recent ones; therefore most of the glacial features currently left on the foreland present a record of the last major glaciation of the Alps, dating from the Last Glacial Maximum (LGM).The timing of the LGM in the Alps (Ivy-Ochs et al., 2008;Monegato et al., 2007) is in good agreement with the maximum expansion of continental ice sheets recorded by the Marine Oxygen Isotope Stage 2 between 29 000 and 14 000 years before present (BP) (Lisiecki and Raymo, 2005).
Based on the interpretation of moraines, trimlines and erratic boulders, Jäckli (1970) first compiled the geomorphological information available on the length and height of former glacial lobes to create a map showing the LGM ice coverage in Switzerland, which was subsequently improved by Bini et al. (2009).The LGM ice extent in Austria was mapped in detail by van Husen (1987).Several reconstruc-tions from surrounding regions of the Alps, including regional reconstructions of the central Swiss Alps (Florineth and Schlüchter, 1998), the Grisons (Florineth, 1998), the Upper Valais to Savoy (Kelly et al., 2004), the French Alps (Buoncristiani and Campy, 2011), the Po plain in Italy (Castiglioni et al., 1999) and Slovenia (Bavec et al., 2004), complete the picture of LGM glaciation extent.They have been compiled in a map covering the entire Alps (Geologische Bundesanstalt, 2013).
Climate information in the LGM was obtained from several studies based on geomorphological, biological information as well as on climate and mass balance modelling.The present-day precipitation pattern in the Alps is dominated by a westerly to northwesterly regime along the northern and western slope (more than 50 % of all days) and by a Mediterranean regime (advection from the southwestern sector) in the south (25 % of all days).However, several sources support an LGM precipitation pattern dominated by southwesterly advection of atmospheric moisture.Based on geomorphological evidence, Florineth and Schlüchter (1998) have shown that the LGM alpine ice cap consisted of at least three major ice domes, all located south of the principal alpine weather divide.The location of the ice domes hints at an LGM precipitation regime more strongly dominated by southwesterly advection than today.An oxygen iso-P.Becker et al.: Numerical reconstruction of the precipitation pattern in the Alps at the Last Glacial Maximum tope record in speleothems from the "Sieben Hengste" cave system in the Bernese Alps (Luetscher et al., 2015) also indicates a preferential southwesterly advection of moisture during the LGM (26 500-23 500 years BP).Furthermore, estimates of basal shear stresses and mass balance gradients by Haeberli and Penz (1985) suggest a strong gradient of climate conditions between the northern and the southern Alps.
A palaeoclimate reconstruction based on pollen records from non-glaciated regions resulted in an annual temperature reduction of 12 ± 3 • C and a reduction in precipitation by 60 ± 20 % north of the Mediterranean Sea (Peyron et al., 1998).Further, a study based on a fossil coleoptera assemblage from Marine Oxygen Isotope Stage 3 by Jost-Stauffer et al. (2005), resulting in mean July temperatures of 12 to 13 • C and mean January/February temperatures of −15 to −7 • C, illustrates the high seasonality of temperatures during the last glacial cycle.
Temperature reductions inferred from atmospheric circulation models are typically less dramatic than reconstructions from proxy records (Kageyama et al., 2001).However, much of this discrepancy has been resolved by accounting for the effect of decreased CO 2 concentrations in the LGM on plant growth through inverse vegetation modelling (Ramstein et al., 2007;Wu et al., 2007).
Using a regional climate model, Strandberg et al. (2011) performed numerical simulations of the climate at LGM.The model results include annual mean air temperatures of 6 to 9 • C below present-day temperatures for the ice-free areas of central Europe, increased winter precipitation rates south and decreased rates north of the Alpine range.However, the model resolution of 50 km is not sufficient to resolve local details of precipitation change in inner-alpine areas.The increase in winter precipitation over southern Europe during the LGM can be explained by a southward displacement of the North Atlantic storm track caused by the high surface elevation of the Laurentide Ice Sheet (Hofer et al., 2012).Indeed, the presence of the Laurentide Ice Sheet induces an enhancement of stability of the North Atlantic eddy-driven jet in this southward shifted position (Merz et al., 2015).
A climate reconstruction of the central European uplands including the Black Forest and the Vosges Mountains based on glacier mass balance modelling (Heyman et al., 2013) compares the modelled area of positive surface mass balance with geomorphological evidence for LGM glaciers.The study concludes that during the LGM, the precipitation was reduced by 25 to 75 %, and the temperature was depressed by 8 to 15 • C compared to today's conditions.Allen et al. (2008a, b, c) applied a glacier-climate model based on mass balance modelling to reconstruct the spatial distribution of glaciated area of the Alps in the LGM to derive a cooling of between 12 and 17 • C north and south of the Alps, respectively.
However, these studies do not include a glaciological model to account for the ice dynamics.Here we use a numerical glacier flow model constrained by the geomorphological ice extent to simulate the entire alpine ice cap under two distinct precipitation patterns.Our specific aim is to analyse the compatibility of the mapped LGM ice cap margin with the present precipitation pattern over the Alps and a modified precipitation pattern to account for prevailing southwesterly moisture advection.
As a numerical model we use the Parallel Ice Sheet Model (PISM, the PISM authors, 2015) to simulate the LGM extent of the alpine ice cap under a constant climate (temperature and precipitation).First, we apply different temperature cooling and precipitation pattern correction factors relative to the present-day climate provided by WorldClim (Hijmans et al., 2005).Second, we introduce a different precipitation pattern, including a reduced precipitation rate in the northern Alps relative to the southern Alps.

Alpine orography and LGM ice extent
Located in central Europe, the Alps form an elongated arc, covering the region between 44 and 47 • N latitude and between 4 and 13 • E longitude.In the north they are bordered by the Jura, Vosges, Black Forest and Bohemian Forest mountain ranges.In the south they arc around the Po plain and the Mediterranean Sea in the southwest.Generally, the Alps are subclassified into the eastern and western Alps.The highest mountains and deepest valleys are located in the western Alps.More precisely, the highest summits are Mont Blanc (4810 m a.s.l.) and Monte Rosa (4634 m a.s.l.).Both are located in the western Alps.In the eastern Alps, Piz Bernina (4049 m a.s.l.) and Ortler (3905 m a.s.l.) are the highest summits.Most valleys are oriented parallel to the alpine arc and slopes are typically steeper in the southern part of the Alps than in the northern part.
Typically, the western Alps receive more precipitation because of an increasing continentality towards the east.During the last glacial cycle, broad areas of the Alps were covered by a contiguous ice cap (Fig. 1, Geologische Bundesanstalt, 2013).This ice cap extended from the Jura Mountains in the west to the Klagenfurt Basin, close to the Slovenian border in the east, and from the area of the present-day Lake Constance in the north to Lake Garda in the south.Smaller ice caps were found on the easternmost reliefs; the largest was located on Hochschwab (van Husen, 2004).
In this study we focus on four of the well-developed major glacier lobes selected in order to represent the ice cap extent in each geographic direction (Fig. 1).The maximum extent of these lobes, which is well constrained by geomorphological evidence, was used to calibrate the model (Sect.3.2).
In the north, our focus is on the Rhine Glacier.It was fed mainly by the Vorderrhein ice dome but also by the Engiadina ice dome (Florineth, 1998) and smaller tributary glaciers draining the northern part of the Alps, and reached a maximum extent near the Danube drainage divide, near Lake Constance.Hence, the Rhine Glacier is mainly affected by the climate north of the Alps.In a geomorphological reconstruction of the chronology of the Rhine Glacier, Keller and Krayss (2005) dated its advance from the confluence of Vorderrhein and Hinterrhein to the LGM extent from 30 000 to 25 000 yr BP.In the east, the focus is on the glacier outflow of the Salzach Glacier.This is the easternmost glacier lobe that advanced into the northern alpine foreland in the LGM.It started retreating at about 21 000 yr BP (van Husen, 1997).In the south, the focus is on the Garda glacier lobe.The Garda Basin acted as a catchment for several glaciers draining the southern part of the Alps such as the Rhaetian Glacier.Several Italian glacier lobes reached their most advanced position of the last glacial cycle between 27 000 and 21 000 yr BP (Ivy-Ochs, 2015).In the west, we study the glacier lobe of Lyon.Several smaller glaciers from the Mont Blanc area such as the Isère and the Arve glaciers contributed to the extent of this lobe (Coutterand et al., 2009).It reached its maximum extent prior to 26 000 yr BP (Buoncristiani and Campy, 2004;Monjuvent and Nicoud, 1988).

Methodology
For numerical simulations of the glaciation of the Alps in the LGM we used PISM (Version 0.6.2,see the PISM authors, 2015).PISM combines submodels to compute ice dynamics and temperature (Sect.3.1), the surface mass balance (Sect.3.2) based on the climate forcing (Sect.3.3) and the lithospheric response.The latter was computed by PISM's Earth deformation model following Lingle and Clark (1985) and Bueler et al. (2007).Our simulations were carried out at a horizontal resolution of 5 km using 180 × 120 grid points.

Ice thermodynamic model
The Stokes model (Greve and Blatter, 2009, chap. 5) is the most physically accurate model for simulating ice dynamics.However, since numerical implementations of this model are computationally very demanding, we used PISM, a simplified hybrid-type model that combines the Shallow Ice Approximation (SIA, Hutter, 1983) and Shallow Shelf Approximation (SSA, Morland, 1987) as a "sliding law".This hybrid scheme is computationally less expensive than the Stokes equations, making it appropriate for the long-term and largescale simulations presented here.Golledge et al. (2012) have successfully used PISM to model the LGM glaciation in the New Zealand Southern Alps, which have a topography similar to that of the European Alps.
The ice velocity was derived by superposing SIA and SSA velocities (Winkelmann et al., 2011).The model also computed the englacial temperature field using an enthalpy-based scheme (Aschwanden et al., 2012) that is capable of conserving energy.Internal energy is given by the temperature of cold ice and the latent heat of the liquid water stored in temperate ice.External energy is provided by the geothermal heat flux at the lowest of 200 lithospheric layers, assumed to be 75 W m −2 (according to Shapiro and Ritzwoller, 2004), P. Becker et al.: Numerical reconstruction of the precipitation pattern in the Alps at the Last Glacial Maximum and constrained by the air temperature from the climate forcing at the surface.
To model basal sliding we used a pseudo-plastic sliding law implemented in PISM (Schoof, 2006).The yield stress was determined by applying a Mohr-Coulomb criterion using an assumed till friction angle of φ = 45 • .Effective pressure, N till , accounted for the till deconsolidation model following Tulaczyk et al. (2000).Additional details on model physics and parameter values not given here can be found in Seguinot et al. (2016, Sect. 2 and Table 1).

Surface mass balance model
To simulate surface ablation and accumulation processes, we used the positive degree-day (PDD) model implemented in PISM, which is an empirical relationship between surface air temperature, precipitation and melt rate (Hock, 2003).
The ablation was assumed to be proportional to the sum of positive temperature degree days.To simulate daily variations in temperature, the PDD model added a variability following Calov and Greve (2005).This temperature variability was assumed to be distributed normally around zero and expressed by a given standard deviation, σ (Calov and Greve, 2005).Since the simulated surface mass balance is sensitive to spatial and seasonal variations in σ , we used σ values computed from the ERA-Interim reanalysis (Dee et al., 2011;Seguinot, 2013).The PDD factors are assumed to be 3 mm • C −1 day −1 for snow and 4 mm • C −1 day −1 for ice.The values are within the range of 3.5 ± 1.4 to 4.6 ± 1.4 mm • C −1 day −1 obtained by Braithwaite (2008), analysing the accumulation-temperature relation with a PDD model over the Greenland ice sheet.The influence of the PDD factors will be analysed in detail in Sect.5.2.
The accumulation rate was given by the solid precipitation rate corresponding to the prevailing precipitation when precipitation was interpreted as rain and discarded; therefore no accumulation occurred.Between T = 0 • C and 2 • C, a linear transition was applied.Possible refreezing of rainwater and melted ice and snow was not taken into account.Finally, the surface mass balance was provided by the difference between surface accumulation and surface ablation.

Climate forcing
To simulate the prevailing climate conditions in the LGM, we applied temperature and precipitation corrections to presentday climate data provided by WorldClim (Hijmans et al., 2005).More precisely, WorldClim provides monthly distributions of temperature and precipitation, which have been generated by a spatial interpolation between weather stations and averaged for the 1950-2000 period.
In a first set of simulations we lowered the temperature field by temperate offsets, following the previous reconstructions from mass balance modelling (Heyman et al., 2013) and pollen records (Peyron et al., 1998).To reduce the precipitation rate we introduced a correction factor c ∈ {0.2, 0.4, 0.6}.Hence, the precipitation P (x, y) at point (x, y) was given as We remark that for all simulations shown in this work, a vertical temperature lapse rate gradient of γ = 6 • C km −1 has been introduced to account for deviations in temperature due to the elevation.
In a second set of simulations, we made two assumptions.First, we defined a weather divide representing a topographic barrier to separate the domain into two areas.We defined a northern domain, N , that received precipitation advected mainly by northwesterly winds from the Atlantic Ocean, and a southern domain, S , that received precipitation advected mainly by southwesterly winds from the Mediterranean Sea.To draw a border between these two domains, we exploited a well-observed weather phenomenon known as "foehn" that divides the central Alps into a northern and southern domain.The (southern) foehn situation occurs when winds come from the south, and it is characterized by moderate to heavy precipitation in the southern Alps and nearly dry, strong foehn winds north of the Alpine Divide (Steinacker, 2006).We assumed that the boundary of the weather divide over the Gotthard and Brenner passes is as shown in Fig. 1 (dashed line).The introduction of the weather divide enabled us to prescribe different precipitation rates in the northern and southern domains.
We assumed a precipitation rate P (x, y) at a point (x, y) of where c N and c S denote correction factors for the northern domain, given by N , and the southern domain, S .The correction factors were given as Here, we assumed a weighting of the precipitation rate in the southern domain S by a factor of 7/5 and in the northern domain N where it was reduced by a factor of 3/5.The values of these weighting factors were obtained by test simulation runs (not shown here) which yielded reasonable results.As a second assumption, we introduced a drying correction factor, to mimic drying of the atmosphere for a given temperature cooling T ∈ {−10, −12, −14} • C. Since this factor is unknown, it will be varied by d ∈ {0.2, 0.4, 0.6}.These values are on the order of 0.5, proposed by Ohmura et al. (1992), which means a halving of the precipitation by decreasing the temperature by 10 • C. The introduction of this drying factor allowed the comparison of the resulting climate forcing with temperature and precipitation values obtained by previous studies (see Sect. 5) due to an expected drying of the atmosphere caused by lowering the temperature of the atmosphere.Figure 2 illustrates the applied precipitation patterns for both the present-day and the modified precipitation pattern.

Model initialization and time evolution
As an initial basal topography, we used the present-day surface elevation given by the SRTM data provided by World-Clim (Hijmans et al., 2005).Our computations were carried out using a horizontal resolution of 5 km.The impact of the resolution is analysed in Sect.5.6.Starting from ice-free conditions, we applied a constant climate forcing over the entire simulation run.Simulations were stopped for analysis as soon as the total ice volume was stabilized.The four reference glaciers often reached a maximum extent several times before the simulations were stopped.However, the time span required to reach a stable total ice volume was highly dependent on the prescribed precipitation and air temperature forcing.It occured between 8000 and 24 000 model years, depending on the specific climate scenario.

Evaluation of steady and transient states
For the evaluation of model results we considered two cases, the steady states, as defined in Sect.3.4 and the transient states corresponding to the ice build-up phase.In each case, we evaluated the performance of the model results by comparing the modelled spatial extent of the ice cap with the ge- omorphologically based map from the Geologische Bundesanstalt (2013).
In the first case, we compared the steady-state extent with the geomorphological map.In particular, we compared the modelled glacier lobe extent of Rhine, Salzach, Garda and the lobe of Lyon (see Fig. 1) to the extents provided by the map.This was done for both the present-day and modified precipitation patterns presented in Sect.3.3.
In the second case, we evaluated the transient deviation of the lobe lengths from the geomorphological ice extent given in the map from the Geologische Bundesanstalt (2013) during the ice build-up.This deviation was measured along the paths of the valleys corresponding to the different lobes.
In this section, we first assume the present-day precipitation pattern to study the sensitivity of glaciation to different climate scenarios (Sect.4.1).Next, we present results obtained using different precipitation correction factors north and south of the Alps (Sect.4.2).

Using present-day precipitation patterns
In order to determine the sensitivity of the modelled glacial extent, we conducted an ensemble of nine simulations using different climate parametrizations.For each cooling T ∈ {−10, −12, −14} • C we considered reduced precipitation rates, given by correction factors of c ∈ {0.2, 0.4, 0.6}.

Steady states
Figure 3 shows the distribution of ice as it reaches a steady state (Sect.3.5).The modelled ice cap is highly sensitive to the applied cooling and precipitation factors.As expected, we observed an increase in glacier extent with decreasing temperature and increasing precipitation rate.
As input climate conditions become more favourable to glaciation, modelled ice coverage appears first in the high mountains in the central part of the Alps (Fig. 3a), then on the northern alpine foreland (see Fig. 3e), and finally, once the northern part of the modelling domain is completely covered, on the southern alpine foreland (Fig. 3i).Using a uniform correction of precipitation, we do not find a climate parametrization that leads to a modelled glacier extent in agreement with the geomorphological reconstruction both north and south of the Alps: in all nine simulations, the geographic centre of the ice cap is modelled further north than in the geomorphological reconstructions.
Table 2 underlines these observations quantitatively, listing the differences in length between the computed and geomorphologically reconstructed termini of the Rhine Glacier (δ R ), the Salzach Glacier (δ S ), the Garda Glacier (δ G ) and the lobe of Lyon (δ L ) as the steady state is reached.A comparison of the length deviation of these individual lobes reveals a different model response between the Rhine and Salzach glacier lobes in the northern part of the model domain, and the Lyon and Garda glacier lobes in the southern part of the model domain.For instance, this different model response can be seen for the median scenario using T = −12 • C, c = 0.4 (Fig. 3e).For this scenario, the Rhine and Salzach glaciers extend far north beyond the geomorphologically mapped ice extent, while the Garda and Lyon glaciers do not reach it.This different model response between the Rhine/Salzach and the Garda/Lyon glaciers occurs for several other climate scenarios as well.
For each scenario used, the root mean square error (RMSE), i δ 2 i , of the differences in length between the modelled and geomorphologically reconstructed glaciers gives a measure of the misfit.The RMSE is minimal for a scenario using a cooling of T = −10 • C and a precipitation correction of c = 0.4 (Fig. 3b; Table 2, second row).However, the deviation in length from the geomorphological reconstruction of the Garda and the Lyon lobes remains significant.
To quantify the misfit between the northern and eastern modelled ice cap margins and the geomorphologically reconstructed ones, we computed northward and eastward shifts as follows.The northward shift (Table 2, N-shift) was computed as the difference, δ R − δ G , between the deviation of the simulated extent of the Rhine and Garda lobes to the geomorphologically based extent.Similarly, the eastward shift (Table 2, E-shift) was computed as the difference in deviation between the Salzach and Lyon lobes, δ S − δ L .While this simple method obviously cannot cover the complexity of the modelled ice cap, it is at least able to give a relative measure.
In all nine climate scenarios considered, we obtained a mean northward shift of more than 64 km and a mean eastward shift of more than 113 km (Table 2, last row).Thus, on average, our simulations produced an ice margin that is shifted many kilometres to the northeast as compared to the geomorphological reconstructions.

Transient states
To clarify the question as to whether the present-day precipitation pattern can explain the geomorphologically obtained LGM extent prior to the modelled steady state, we evaluated the transient evolution of the four reference glacier lobe deviations relative to the geomorphological LGM extent, δ i (t) (Fig. 4).In two scenarios, using T = −12 • C and c = 0.4 and c = 0.6 (Fig. 3e, d), as well as the three scenarios using T = −14 • C (Fig. 3g, h, i), the modelled ice cap mostly reaches the LGM extent.For these climate scenarios, a distinct minimum of the RMSE function given by i δ 2 i can be identified.The corresponding time series of lobe deviations, δ R (t), δ S (t), δ G (t), δ L (t) (Fig. 4), indicate a constant delay between the lobes of Rhine and Salzach, on the one hand, and the lobes of Garda and Lyon on the other hand.Hence, a good fit between modelled ice margin and the geomorphological LGM extent is not reached simultaneously for all lobes in any climate scenario.

Using modified LGM precipitation patterns
To improve the partitioning of the modelled ice extent between the northern and the southern parts of the model domain, we applied a modified precipitation pattern, weighting the precipitation rate in the southern part of the Alps, as presented in Sect.3.3.Table 2. Deviation of the modelled length of selected lobes, δ i , relative to their geomorphological reconstruction, using the same precipitation correction factor c, north and south of the Alps.The northward shift of the ice cap (N-shift) was computed as the difference, δ R −δ G , between the Rhine and Garda lobe length deviations.The eastward shift of the ice cap (E-shift) was computed as the difference, δ S − δ L , between the Salzach and Lyon lobe length deviations.Bold text indicates the scenario with minimal root mean square error (RMSE), given by i δ 2 i .The median scenario using T = −12 • C and c = 0.4 is highlighted in italic.

Steady state
Figure 5 illustrates the modelled glacier extent at steady state, obtained using the corrected precipitation pattern.Temperature was again varied by T ∈ {−10, −12, −14} • C. The precipitation was reduced by weighting factors c N and c S as described in Sect.3.3.The modelled glaciated area (Fig. 5) is now generally displaced southwards as compared to presentday precipitation pattern scenarios (Fig. 3), resulting in a smaller misfit to the geomorphological reconstruction.A more balanced glacier expansion in the north and the south was obtained for most scenarios.
The deviation of the modelled glacier lobes relative to the geomorphologically reconstructed glacier termini is presented in Table 3.As a matter of fact the misfit in glacier expansion between the northern and the southern parts of the Alps was improved (Table 3, N-shift, E-shift) by using the modified precipitation pattern.This is shown by the average northward shift that was diminished from 64 km (Table 2, last row) to 6 km (Table 3, last row).Moreover, the average eastward shift is diminished from 113 km (Table 2, last row) to 76 km (Table 3, last row).
We obtained the best agreement between the modelled lobe extents and the geomorphological reconstructions for the climate scenario using T = −14 • C, c N = 0.17 and c S = 0.39.However, this scenario did not resolve the shape of single lobes with the same quality as the median scenario of T = −12 • C, c N = 0.20 and c S = 0.47.In general, however, we find a striking improvement in the fitting of both in the north and the south using modified precipitation patterns compared to the present-day precipitation patterns.

Transient states
Figure 6 illustrates the transient evolution glacier lobe length deviations using the modified precipitation pattern.In contrast to climate scenarios using the present-day precipitation pattern (Fig. 4), the delay between the Rhine and Salzach lobes on the one hand and the Garda and Lyon glacier lobes on the other hand was significantly decreased.Furthermore, the minimum values of the RMSE functions were smaller.Hence, a better fit between modelled ice margin and the geomorphological LGM extent was reached more simultaneously for all lobes when using the improved precipitation pattern weighting the precipitation rate in the southern Alps.

LGM climate in the Alps
Simulations using the modified precipitation pattern result in a better fit between the modelled LGM alpine ice cap and the geomorphological reconstruction.This modified precipitation pattern contains two corrections.First, we reduced the precipitation rate according to the decrease in air temperature to simulate an atmospheric drying effect.Second, we redistributed the precipitation from the northern to the southern part of the Alps to reflect a southward shift of the North Atlantic storm track.The results obtained from our simulations support the geomorphological indications discussed by Florineth and Schlüchter (1998) of an LGM precipitation pattern that was dominated by south-southwestern winds, which is in contrast to the present-day precipitation pattern that is dominated by westerly to northwesterly winds.
Our best-fitting climate scenario is given by a cooling of T = −12 • C and a precipitation rate of 20 % in the northern Alps and 47 % in the southern Alps as compared to that of today.These values are comparable to those obtained by Heyman et al. (2013), who derived a precipitation reduction of 25-75 % and a cooling of 8-15 • C for small glaciers in central Europe, north of the Alps.A study of pollen records by Peyron et al. (1998) indicates a temperature cooling of 12 ± 3 • C which are similar to our results.Atmospheric models indicate a less severe cooling (Ramstein et al., 2007).For instance, Strandberg et al. (2011) model temperatures of 6-9 • C below present-day values over the ice-free areas which is above our modelled temperature.However, this study indicates a south-dominated precipitation pattern, which is similar to our simulations.

Sensitivity to PDD factors
PDD factors are dependent on the surface energy balance, which is in turn governed by environmental factors such as albedo, wind speed and temperature variability (Lang and Braun, 1990;Ohmura, 2001).These variables are highly uncertain in the course of the last glaciation.Albedo, for instance, is highly dependent on the dust deposited on the glaciers, which probably varied according to the climate over the time during the last glacial cycle.A study by van den Broeke et al. (2008) presenting field observation data from weather stations in Greenland indicates a wide range of possible PDD factors during the LGM.Since PDD factors during the LGM are uncertain, additional simulation experiments (Fig. 7) demonstrate the sensitivity of our model results with respect to the PDD factors.Here, the modified LGM precipitation pattern, described in Sect.3.3, is used.The change in the PDD factors by ±1 mm • C −1 day −1 corresponds to a change in cooling by ±0.5 to 1.5 • C.
The relative distribution of ice during the ice build-up phase is independent of the PDD factors.Only the total amount of ice is affected by changing these parameters.Thus we cannot conclude that the best matching climate scenario using a cooling of T = −12 • C determines a realistic value for the temperature in the LGM.However, the spatial distribution of ice on a large scale is not affected by the choice of PDD factors.Hence, the choice of PDD factors does not influence the precipitation pattern itself.

Remaining discrepancies of the modelled ice extent
Our preferred climate scenario, the median scenario (Fig. 5e), uses T = −12 • C, c N = 0.20 and c S = 0.47.Although the scenario using T = −14 • C, c N = 0.17 and c S = 0.39 (Fig. 5h) results in a lower RMSE between the modelled and geomorphological reconstructed lobes of Rhine, Salzach, Garda and Lyon (Table 3), our preferred scenario shows a better resolved shape of the alpine ice cap.It is in good agreement with the geomorphological reconstruction along the southern margin of the modelled alpine ice cap.However, significant discrepancies between the model results and the mapped ice cap margins remain in other geographic directions.
Along the northern margin, the alpine topography is given by a flatter relief compared to the south.In this case, a change in the climate conditions strongly affects the length of the north-oriented glacier lobes.In the northern alpine foreland, the length of the glacier lobes between Rhone and Salzach agrees well with the mapped lobes.However, we identify a significant discrepancy of the ice coverage of the Napf region showing a strong tendency to glaciate, which is in contrast to the geomorphological reconstruction (Bini et al., 2009).Significant disagreement of the modelled LGM extent to the mapped reconstruction is also obtained in the east, especially in the region of the Drau Glacier in the Klagenfurt Basin and the Hochschwab Glacier.Both the model and the smoothed margin of the mapped LGM extent are not capable of resolving the patchy ice coverage, mapped in detail by van Husen (1987).Further, the relief in the east is flatter than in the west, www.geogr-helv.net/71/173/2016/making the glacier extent more sensitive to climate conditions as shown by single panels in Fig. 5.

Rhine Salzach
In contrast to glacier mass balance modelling (Heyman et al., 2013), the present results indicate a strong glaciation in the east.Hence our results do not support the theory of a temperature anomaly gradient increasing from east to west as discussed by Heyman et al. (2013).However, it is possible that the choice of the border of the assumed weather divide has hidden this effect.
Overall, in the south, the steep relief leads to robust matching of the LGM ice extent pattern.Nevertheless, certain scenarios tend to show a large glacier lobe advancing from Monte Rosa to the Piedmont region.
In the southwest we observe a more pronounced disagreement of the glacier extent in the Maritime Alps, calling into question the simple approach that consists of only two weather zones for this area.The glaciation of the Jura Mountains is in good agreement with the geomorphological evidence.We note that the reason that the Vosges and Black Forest mountains are not covered by ice is the smoothing of the topography due to an insufficient grid resolution.

Timing of glaciation
For the purpose of comparing the simulated glacier extents with the geomorphological reconstruction, it is essential to pay attention to the growth of the ice caps in nature and in the simulation.In nature the process of glaciation during the LGM is a transient process with a changing glaciated area, albedo, energy balance, surface elevation and sedimentation, www.geogr-helv.net/71/173/2016/Geogr.Helv., 71, 173-187, 2016 etc.By contrast, our simulated ice caps represent a product forced by a constant climate to reach a steady state between 8000 and 24 000 model years.Hence the inferred climate conditions can only be understood as a mean of the timespan between the final onset of cooling and the LGM.During this period, the glaciers experienced multiple fluctuations in extent due to the volatile climate conditions during the last glacial cycle with an ultimate decay of the glacier lobes beginning before about 20 000 years BP.
Beside steady states, we also investigated the temporal evolution of the ice build-up process to clarify the question of the possible existence of transient states that match the mapped LGM extent using present-day precipitation.Fig. 4 illustrates the deviations from the mapped LGM extent of the reference glaciers which are not vanishing isochronologically for any modelled climate scenario.Hence we cannot identify a suitable transient state, matching the LGM extent while using the present-day precipitation pattern.By contrast, there is a constant delay between the advancing northern-affected glaciers of Rhine and Salzach and the delayed southern-affected glaciers of Garda and Lyon throughout all climate scenarios.This delay of lobe length deviations between the northernand southern-affected glaciers disappears using the improved precipitation pattern (Fig. 6).Hence transient states matching the LGM extent of the reference glaciers can be found using the improved precipitation pattern (e.g. for the climate scenario using T = −14 • C, c N = 0.29 and c S = 0.69 after approximately 3500 model years; see Fig. 6i).Furthermore, we compare the speed of the advancing Rhine lobe with the geomorphological evidence.Figure 6 illustrates a duration of less than 5000 years for the maximum advance of the Rhine Glacier (at the Danube drainage divide near Randen) for all climate scenarios in which the mapped LGM extent is reached (Fig. 6e, f, h, i).This agrees well with the geomorphological chronology presented by Keller and Krayss (2005) (see Fig. 9a), who have estimated a duration of the lobe advance between 4000 and 5000 years.

Ice thickness
Our simulations indicate a systematic overestimation of ice thickness compared to the map of Bini et al. (2009) for all climate scenarios.The most important inconsistencies were found in the accumulation area, especially the upper Rhone valley and the Vorderrhein valley, which typically overestimate it by more than 500 m.This overestimation in ice thickness is caused by the obtained ice-cap-like glacier surface geometry.It is obtained for all tested climate scenarios as well as for all test cases, using different sliding conditions (not presented here) or climate conditions.The maps by Bini et al. (2009) or Jäckli (1970) show a more network-like glacier flow that is governed by the topography of the underlying valleys.
First, it might be possible to explain this discrepancy by a change in climate conditions at the time of ice cap formation.In nature, an ice cap could potentially induce a different climate in its interior characterized by low precipitation rates and strong winds redistributing the snow, causing a local decrease in the accumulation rate at high elevation.By contrast, our model is forced over the entire simulation run by a constant accumulation rate and does not account for snow redistribution processes.
Second, the sedimentary infill in the valleys before the LGM onset is unknown.Geological studies indicate a sedimentation thickness of more than 1000 m in the Rhone valley (see Stucki and Schlunegger, 2013, Fig. 4) in the present day.The height of the infill could also affect the ice thickness.However, the sedimentation infill at LGM is unknown and test cases did not show a significant sensitivity of the ice thickness to an infill-free bedrock topography.
Third, the geomorphological mapping information of ice thickness in the LGM is based on glacial trimlines.However, trimlines can only give a minimum value for the ice thickness in a given area.

Sensitivity to horizontal resolution
In the model the computational grid size is a source of systematical error, thus limiting the validity of the results.The 5 km resolution of the digital terrain model used is capable of resolving the main characteristics of the alpine landscape such as large valleys and mountain ranges that dominated the ice flow.Test cases using a 2 km resolution (not presented here) showed slightly extended lobes (< 5 % of total length) and a more distinct shape of the lobes due to the resolution.This means the extents of the obtained glacier lobes computed on a 5 km grid are typically slightly underestimated in length and slightly overestimated in width.

Conclusions
In this study, simulations of the LGM alpine ice cap were run using a numerical ice-flow model forced by present-day and modified precipitation patterns.Simulations forced by the present-day precipitation pattern systematically resulted in an extensive glaciation in the northern alpine foreland relative to glacier extent in the southern Alps.In fact, only the introduction of decreased precipitation in the northern Alps yielded a modelled ice cap extent comparable to the geomorphological reconstructions.Despite remaining discrepancies in modelled ice extent and thickness, this result demonstrates the compatibility of the mapped LGM ice cap margin with a different atmospheric circulation pattern than the present-day one.Our simulations support an LGM dominated by southwesterly moisture advection to the Alps, and thus a southward shift of the North Atlantic storm track during the LGM.The present results have been restricted to the specific limitations of the model.These limitations relate to the model resolution and parametrizations such as ice softness and sliding parameters.For future studies we recommend that they be further investigated and improved upon.In particular, the way in which the ice thickness deviates from the geomorphological evidence should be subjected to detailed scrutiny.

Data availability
Data are not publicly accessible.

Figure 1 .
Figure 1.Simplified LGM glacier extent in the Alps, redrawn based on the map from the Geologische Bundesanstalt (2013).Arrows indicate the flow direction of the four glaciers used to quantify the ice cap extent in each geographic direction.(MB indicates the Mont Blanc range and MR the Monte Rosa range, showing the two main accumulation areas in the western Alps.PB indicates the Piz Bernina and O the Ortler range, showing the two main accumulation areas in the eastern Alps.The dashed red line marks the assumed weather divide.)

Figure 2 .
Figure 2. Applied mean monthly precipitation rate patterns for January and July in m yr −1 for the unimproved, present-day precipitation pattern and the modified precipitation pattern (second row).The red line indicates the assumed weather divide that is used in the modified precipitation pattern to separate the domain into two areas.

Figure 3 .
Figure 3. Modelled extent of the LGM alpine ice cap under different temperature cooling, T , and different precipitation correction factors, c.The red line indicates the geomorphologically based glaciation extent (Fig. 1).Contour lines correspond to 3000, 2000 and 1000 m a.s.l.

Figure 4 .
Figure 4. Time series of the deviation between modelled lobe length of selected glacier lobes (Sect.3) and the geomorphological LGM extent, δ i (t), (solid colour lines, left axis) in the case of forcing the model with present-day precipitation patterns.A negative deviation indicates a simulated lobe shorter than the field-based measurement.A positive deviation indicates a glacier lobe extent exceeding the fieldbased measurement.The root mean square error, RMSE, i δ 2 i of the deviations (dashed black line, right axis), represents a measure of the misfit between the simulated glacier lobes and the field-based evidence of the LGM extent.

Figure 5 .
Figure 5. Modelled extent of the LGM alpine ice cap under varied temperature coolings, T and different precipitation rates north and south of the Alps, c N and c S .The red line indicates the geomorphologically based glaciation extent (Fig. 1).Contour lines correspond to 3000, 2000 and 1000 m a.s.l.

Figure 6 .
Figure 6.Time series of the deviation between modelled lobe length of selected glacier lobes (Sect.3) and the geomorphological LGM extent, δ i (t), (solid colour lines, left axis) in the case of forcing the model with modified LGM precipitation patterns (Sect.3.3).A negative deviation indicates a simulated lobe shorter than the field-based measurement.A positive deviation indicates a glacier lobe extent exceeding the field-based measurement.The root mean square error, RMSE, i δ 2 i of the deviations (dashed black line, right axis), represents a measure of the misfit between the simulated glacier lobes and the field-based evidence of the LGM extent.

Figure 7 .
Figure 7. Averaged precipitation factors, (c N + c S )/2, corresponding to a minimal misfit between the modelled steady-state extent of the Rhine, Salzach, Garda and Lyon lobes and their geomorphological reconstruction.The dark grey line corresponds to our default choice of PDD factors of 3 mm • C −1 day −1 for snow, and 4 mm • C −1 day −1 for ice.The other lines correspond to PDD factors increased and by 1 mm • C −1 day −1 relative to these default values.

Table 1 .
Values and symbols of the physical quantities used in the model.DDF denotes degree-day factor.