The glacier model
The evolution of glacier ice thickness, denoted as \(h(x,y,t)\), starting from an initial glacier shape, is governed by mass conservation, which connects elevation change, ice dynamics, and mass balance through the following continuity equation:
$$\frac{\partial h}{\partial t}+\nabla \cdot \left(\bar{{{{\bf{u}}}}}h\right)={{{\rm{MB}}}},$$
(1)
where \(\nabla\) represents the divergence operator for the horizontal variables \((x,y)\), \(\bar{{{{\bf{u}}}}}=\left(\bar{u},\,\bar{v}\right)\) denotes the vertically-averaged horizontal ice velocity field, and \({{{\rm{MB}}}}\) represents the surface and basal mass balance functions. Equation (1) is solved using an explicit upwind finite volume scheme on a regular grid which allows the model to update the ice thickness while conserving mass. In the following sections, we describe in turn individual IGM31,32 sub-models used in this study for simulating processes of ice flow, ice enthalpy, climate forcing, surface mass balance, isostatic adjustment, and avalanching.
Ice flow
Ice flow dynamics is modelled using Glen’s flow law:
where \(D\) and \(\tau\) denote the strain rate and deviatoric stress tensors, with \(A\) representing the rate factor, and \(n=3\) as Glen’s exponent. Here, we use the Blatter-Pattyn model33, which disregards second-order terms in the thickness/length ratio in the momentum conservation equation. This modification makes solving the stress balance easier than with the original Full-Stokes model. For our boundary condition, we use a nonlinear Weertman friction condition (e.g. Schoof & Hewitt73), relating the basal shear stress \({\tau }_{b}\) to the sliding velocity \({{{{\bf{u}}}}}_{b}\) as follows:
$${\tau }_{b}=-c{|{{{{\bf{u}}}}}_{b}|}^{m-1}{{{{\bf{u}}}}}_{b},$$
(3)
where \(m \, > \,0,{c}=c\left(x,y\right) \, > \,0\) is the sliding coefficient.
Rather than using a traditional solver, IGM models the ice flow using a physics-informed deep learning emulator74 trained to minimise the energy associated with the Blatter-Pattyn equation32. Thanks to its efficient evaluation and training on GPUs, the neural network considerably reduces computational cost with negligible accuracy losses relative to a traditional ice-flow solver32. The neural network features 16 two-dimensional convolutional layers representing 140, 000 trainable parameters. We use the hyper-parameters found by Jouvet et al.31. To obtain initial weights and facilitate convergence, the neural network is pretrained over a diverse catalogue of glaciers and flow regimes32. Moreover, it is frequently re-trained during transient IGM simulations (every seven iterations) to adjust to new glacier states obtained through time. This frequency was found to be a good trade-off maintaining accuracy while keeping computational cost relatively low.
Ice enthalpy
Temperature within the ice is modelled with an energy-conservative enthalpy model following Aschwanden et al.34. Ice enthalpy, \(E\), is a function of ice temperature, \(T\), and ice water content, \({{{\rm{\omega }}}}\):
$${E}(T,\omega,p)=\left\{\begin{array}{l}{c}_{i}(T-{T}_{ref}),if\,{T} < {T}_{pmp},\\ {E}_{pmp}+L\omega,if\,T={T}_{pmp},\,0\le \omega \end{array}\right.,$$
(4)
where the temperature, Tpmp, and enthalpy, Epmp, at pressure-melting point of ice are defined by
$${T}_{{pmp}}={T}_{0}-\beta p,$$
(5)
$${E}_{{pmp}}={c}_{i}\left({T}_{{pmp}}\left(p\right)-{T}_{{ref}}\right)$$
(6)
According to the definition of enthalpy prescribed above, we have two possible modes: i) When the ice is cold (i.e. below the melting point), the enthalpy is simply proportional to the temperature minus a reference temperature. ii) When the ice is temperate, the enthalpy continues to increase. In this case, the additional component, Lω, accounts for the creation of water content through energy transfer. The enthalpy model consists of an advection-diffusion equation, with horizontal diffusion being neglected, and strain heating and drainage as source terms. At the modelled ice surface, the enthalpy equation is constrained by the surface temperature provided by the climate forcing. Borehole data show the offset between temperatures of surface air and of the active ice layer at glacier surfaces is challenging to constrain and varies through both space and time75. Here, we set this offset as an ensemble-varying parameter with possible values ranging between 1 °C and 3 °C (surface ice being slightly warmer), bracketing the median (1.85) of observed offset values (n = 41) reported by Zagorodnov et al.75. At the glacier bed, there are multiple boundary conditions for the enthalpy equation depending on the bottom ice layer and bed surface temperature34, the latter being forced, in this study, by geothermal heat flux data from Goutorbe et al.35. The 3D advection-diffusion equation is solved using a semi-implicit scheme with finite differences at each time step, defined by the time stepping of the mass conservation. As a result, the ice enthalpy, temperature, and basal melt rate are all updated at each time step. The enthalpy impacts both internal ice flow by influencing the rate factor, \(A\), via the Glen-Paterson-Budd-Lliboutry-Duval law56:
$$A\left(T,\omega \right)=A\exp (-Q/(R{T}_{{pa}}))(1+181.25\omega ),$$
(7)
and basal sliding, \(c\), via meltwater production when pressure melting point is reached. Our implementation of the enthalpy formulation in IGM successfully passed the two benchmark experiments proposed by Kleiner et al.76 and Hewitt & Schoof77. These tests, as well as more details regarding the enthalpy model can be found alongside IGM’s source code (https://github.com/jouvetg/igm).
Sliding parameterization
Following Bueler & van Pelt78, the basal water thickness in the till layer, \({W}_{{till}}\), is computed from the basal melt rate, \({m}_{b}\), obtained from the enthalpy as follows:
$$\frac{\partial {W}_{{till}}}{\partial z}=\frac{{m}_{b}}{{\rho }_{w}}-{C}_{{dr}},$$
(8)
where \({C}_{{dr}}\) is a simple drainage parameter. The till layer is assumed to be saturated when the basal water layer thickness reaches a caping value of \({W}_{{till}}^{\max }=2{{{\rm{m}}}}\). The effective thickness of water within the till layer \({N}_{{till}}\) is computed from the saturation ratio \(s={W}_{{till}}/{W}_{{till}}^{\max }\) by the formula78:
$${N}_{{till}}=\min \left\{p,{N}_{0}{\left(\frac{\delta P}{{N}_{0}}\right)}^{s}{10}^{({e}_{0}/{C}_{c})(1-s)}\right\},$$
(9)
Where \(p\) is the ice overburden pressure and the remaining parameters are constant. The sliding coefficient, \(c\), in (3) is defined by the Mohr-Coulomb law56 that involves the effective pressure in the till \({N}_{{till}}\):
$$c={\tau }_{c}{u}_{{th}}^{-m}={N}_{{till}}\tan (\phi ){u}_{{th}}^{-m},$$
(10)
where \(\phi\) is the till friction angle25. Here, using the assumption that basal materials are generally weaker (softer sediments) in valley troughs than over mountain tops38, \(\phi\) is parameterised to be a piece-wise linear function of bed elevation, b:
$$\phi \left(x,y\right)=\left\{\begin{array}{c}{\phi }_{\min },\hfill \,b\left(x,y\right)\, \le \, {b}_{\min },\\ {\phi }_{\min }+\left(b\left(x,y\right)-{b}_{\min }\right)M,\quad{b}_{\min }\, < \, b\left(x,y\right)\, < \, {b}_{\max },\,\\ {\phi }_{\max },\hfill \,{b}_{\max }\, \le \, b\left(x,y\right).\end{array}\right.$$
(11)
Between ensemble simulations, we modify the elevation-dependency of \(\phi\) by varying upper and lower bed elevation thresholds (\({b}_{\min }\), \({b}_{\max }\)) between ranges of (−500, −100) metres, and (2400, 3000) metres, respectively, while till friction angle thresholds (\({\phi }_{\min }\), \({\phi }_{\max }\)) are set to values of 15 and 50 (see Supplementary Table 2). To further modify the sliding coefficient for a given basal shear stress, the parameter \({u}_{{th}}\) also varies within our ensemble between 100 and 2000 m yr-1.
Climate forcing
In this study, input climate forcing fields are taken from Jouvet et al.22. First, the Weather Research and Forecasting regional climate model27 was used to downscale time slice simulations of a global Earth system model79 to high-resolution (2 km) over the European Alps. This workflow produced climate snapshots, including weekly mean and standard deviation data, for the pre-industrial (1850 AD), LGM in the Alps (~ 24 ka) and Marine Isotope Stage 4 (65 ka) periods (see Figs. 5–7 in Jouvet et al.22 for more details). Here, we extend climate data between these three snapshot states using a glacial index approach (e.g. Niu et al.41). This method creates continuous climate fields with two given states: the pre-industrial climate snapshot with limited ice cover in the Alps, and the LGM climate snapshot. The glacial index function is built by linearly rescaling a climate proxy signal such that the glacial index is close to 1 at the LGM and close to 0 at the pre-industrial. In this study, we use an Alps-specific climate proxy signal as input to our glacial index scheme (Fig. 2e), which combines the Bergsee lacustrine record42 (35–30 ka: Black Forest, Southern Germany) and the Sieben Hengste speleothem δ18O record43 (30–18 ka: Bernese Alps, Switzerland). Note that input air temperature fields refer to the surface topography given as input to the regional climate model, which features glaciers at their maximum extent for Marine Isotope Stage 4 and the LGM states, and the present-day topography for the pre-industrial state. To simulate the temperature when the modelled surface deviates from the reference one, we apply a vertical and linear correction using an atmospheric lapse rate of 6 °C km−1.
While the above climate forcing improved the overall model-data fit in the LGM extent of the AIF relative to previous studies20,24, certain outlet glacier extents remained either too large or too small despite varying non-climatic parameter extensively. We assume these isolated misfits are likely related to uncertainties in the input climate and/or limitations of the glacial index approach, and thus implement an ensemble-varying and glacier-catchment-specific precipitation offset scheme. To do so, we apply time-independent scalar multipliers to the input precipitation field over the Rhein (reference multiplier = 0.9), Isar (1.33), Jura (0.7), Drau (1.33) and southernmost Alps (0.7) regions. In other regions, precipitation remains equal to the original input field. The magnitude of these offsets is made simulation-dependent by using an ensemble-varying multiplier parameter (ranging between 0.85 and 1.15) that further modifies reference offset values. Original precipitation values can thus be modified by between -40% and +50%, depending on the region and ensemble simulation (See Supplementary Table 1).
Surface mass balance
In this study, we parameterize IGM with a combined snow accumulation and positive degree-day model80 to compute surface mass balance from input temperature and precipitation fields. In this scheme, precipitation generates surface accumulation (falls as snow) when air temperature is below 0 °C and causes no accumulation (falls as rain) when temperature is above 2 °C, with a linear transition in between. Surface ablation, on the other hand, is computed proportionally to the number of positive degree days. The positive degree day integral is numerically approximated using week-long sub-intervals, following Calov & Greve36. Positive degree day parameters are not well constrained and can vary in space and time. Thus, in our ensemble, the melt factor for ice is simulation-dependent and varies between 6 and 9 mm w.e. d−1 °C−1 (see Supplementary Table 1). The melt factor for snow remains constant at 3 mm w.e. d−1 °C−1. Our positive degree day scheme also models the refreezing (turned into net accumulation) of a given proportion of the computed melt. This proportion is here made simulation-dependent and varies between 50% and 70% within the ensemble (see Supplementary Table 1). Note that no proglacial lake module is implemented in this model setup, meaning that all ice is assumed to be land-terminating. We consider this assumption to have little impact on the LGM geometry of the AIF since most overdeepened basins of the Alpine foreland were eventually ice-filled during maximum glacier advance18.
Model initialisation
All IGM ensemble simulations are initialised with ice-free conditions at 35 ka. Although unrealistic, a sensitivity analysis reveals that starting an ice-free, AIF-wide simulation at 40 ka instead does not impact modelling results at the LGM, as diagnostic model variables converge after 4–5 kyr of running the model (see Supplementary Fig. 13). Therefore, as suggested by previous modelling work22, the response time of the AIF to climate change during the last glaciation does not exceed 4–5 millennia.
Model validation
Before applying IGM and conducting AIF-wide simulations at 300 m resolution, we quantitatively assessed the model’s suitability and fidelity in simulating the last glaciation of the European Alps. To do so, we ran simulations at the same 2 km spatial resolution as Jouvet et al.22 to quantitatively compare outputs with PISM25, a widely used and well-tested model25. Unlike our higher-resolution (300 m) ensemble runs (see ‘Results’ section), these test simulations use the same glacial index climate signal (EPICA ice core40), basal topography, and space-independent basal till friction angle (ϕ = 30°) parameterization as Jouvet et al.22. Furthermore, no avalanche module is employed in this 2 km setup. This experiment resulted in a 2 km simulation with IGM producing AIF extents, volumes, ice velocities, surface mass balances and basal conditions that are consistent with PISM outputs from Jouvet et al.22 (see Supplementary Figs. 14–16). Between the two models, maximum AIF volume and areal extent (reached in this case at ~24.5 ka) vary by less than 2%, while LGM ice thickness differences remain minimal (−6 ± 147 m). At 2 km, Alps-wide differences in LGM ice flux between IGM and PISM are minor with surface, basal, and depth-averaged velocities varying by 2 ± 110 m yr−1, -27 ± 110 m yr−1, and −6 ± 108 m yr−1, respectively. This test shows that IGM can be used to model the AIF’s last glaciation with results comparable to a well-established ice-sheet model (see Supplementary Figs. 14–16). Results are not expected to be 100% identical, however, as important differences remain between the two models, the most potent of which is the ice flow stress balance, approximated with the SSA/SIA hybrid model in PISM25, and with the higher-order Blatter-Pattyn33 model in IGM32.
Isostatic adjustment
In the European Alps, the lithospheric deflection caused by ice loading during the LGM was large enough to notably influence ice surface slope and elevation, and thus ice flow dynamics and mass balance (see Supplementary Fig. 12). To account for space- and time-dependent bed deflection, we couple IGM with the gFlex model39 which dynamically computes the flexural isostatic adjustment using the two-dimensional elastic thin-plate Eq. (1) in Mey et al.21. For gFlex-specific calculations, the frequency of iteration is here set to 50 years while the spatial resolution remains at 2 km. At the entire Alps scale, the impact of a higher frequency or spatial resolution on modelled AIF evolution is negligible81. For our 300 m resolution IGM runs, we feed gFlex a space- and time-independent lithospheric effective elastic thickness, representing the resistance to bending under specific vertical loads. This elastic thickness is challenging to constrain and is thus set as an ensemble-varying parameter ranging from 35 to 50 km, after Mey et al.21 (see Supplementary Table 1).
Avalanche scheme
Increasing spatial resolution with IGM produces finer and steeper bed topographies in upper alpine catchments and accumulation zones. Accurate modelling of transient glacier evolution on such topographies requires an approximate representation of avalanching which impacts ice accumulation, surface elevations, and flow velocities. With IGM, we thus make use of an avalanche module that redistributes modelled accumulation downslope until the glacier surface reaches a given angle of repose value, here set to 45° across the domain and for all simulations. A sensitivity analysis on the best-fit simulation (number 37) reveals that using a value of 35°, on the low-end of reported angle of repose values for Alpine glaciers, does not generate a notable impact on AIF-wide model-data fit in LGM ice extent (−0.03%) and thickness (−3%), relative to 45° (see Supplementary Fig. 17). In our simulations, the frequency of the avalanche module updates is set to 5 years.
Empirical trimline elevation data
Field-based estimation of trimline location and elevation is non-trivial, yields geomorphological uncertainties, and can sometimes be challenged by dangerous access conditions forcing remote observations. Therefore, in this study, empirical trimline elevation data (n = 396), gathered from literature10,11,12,13,14,15,16, are quality-controlled by comparing them against independent elevations extracted from high-resolution (≤ 5 m) digital elevation models at their reported locations (see Supplementary Table 3). Given the relatively low (<10 m) vertical error of these satellite-derived products, a trimline data point yielding a >50 m offset between the two independent elevations was here considered an outlier likely due to error in measurement of either trimline geolocation or elevation. Trimlines yielding >50 m offsets were found to represent ~11% of the original dataset (n = 43 out of 396). The impact of removing these 43 presumed outliers on the mean model-data misfit between modelled ice surface elevations and reported trimline elevations is negligible (on average a ~ 2% reduction). However, it notably reduces the standard deviation of this misfit, on average by ~40%, thus decreasing scatter and generating a notable improvement in model-data ice thickness agreement. Details of the different digital elevation models used for this analysis are listed in Supplementary Table 3.
The empirical LGM outline of the AIF
The empirical LGM outline of the AIF used in this study was originally produced by an Alps-wide compilation of geomorphological and geochronological evidence by Ehlers et al.17. Since this key study, the outline has been updated by a series of empirical investigations improving the quality of LGM margin reconstructions in specific sectors of the Alpine foreland. In this work, we use the most up-to-date version of the outline. Here, we provide a summary of the main studies which have updated the LGM AIF outline since Ehlers et al.17. Gianotti et al.50,51 have made updates to LGM margins of the Ivrea outlet glacier in the Dora Baltea region. Braakhekke et al.49 have done so for the Orta region, Kamleitner et al.5 for the Verbano region, Ravazzi et al.53 for the Oglio region, Monegato et al.48 for the Garda region, Federici et al.47 for the Gesso region, Ribolini et al.8 for the Stura region, and Ivy-Ochs et al.52 for the Dora Riparia region.