A data-consistent model of the last glaciation in the Alps achieved with physics-driven AI

Domain: earth_science
Published: 2025-01-20
Editorial outcome: minor_revision

Abstract
25 thousand years ago, the European Alps were covered by the kilometre-thick Alpine Ice Field. Numerical modelling of this glaciation has been challenged by model-data disagreements, including overestimations of ice thickness. We tackle this issue by applying the Instructed Glacier Model, a three-dimensional model enhanced with physics-informed machine learning. This approach allows us to produce 100 Alps-wide and 17 thousand-year-long simulations at 300 m resolution. Previously unfeasible due to computational costs, our experiment both increases model-data agreement in ice extent and reduces the offset in ice thickness by between 200% and 450% relative to previous studies. Our results have implications for better estimating former ice velocities, ice temperature, basal conditions, erosion processes, and paleoclimate in the Alps. This study demonstrates that physics-informed machine learning can help overcome the bottleneck of high-resolution glacier modelling and better test parameterisations, both of which are required to accurately describe complex topographies and ice dynamics.

Introduction
The extent, thickness, and flow geometry of the European Alpine Ice Field (AIF) during Quaternary glaciations has been studied for more than a century1–4. The result is an abundant library of mapped and dated ice-contact glaciogenic landforms and sediments (e.g. moraines and tills)5–8, glaciofluvial deposits9, and trimlines10–16, making the extent of the AIF during the Last Glacial Maximum (LGM: ~30–19 ka17,18) one of the most well known in the world11,18. Yet, numerous controversies remain regarding the spatio-temporal evolution of the AIF’s former thickness, velocity patterns, basal conditions, outlet glacier asynchronies, erosion potential, and isostatic adjustment. Obtaining model simulations of the LGM ice field that match empirical data has implications for understanding such mechanisms and can also improve reconstructions of late-Quaternary human and ecological history19.

Under these motivations, over the past 20 years, several studies20–24 have produced simulations of AIF glaciations using glacier evolution models such as the Parallel Ice Sheet Model (PISM25) and the integrated second-order shallow ice approximation model (iSOSIA26), with the main objective to closely fit moraine and trimline evidence. As a result, model-data misfits in LGM ice extent were incrementally reduced, with the latest study22 making progress by using an improved climate forcing27. However, a major discrepancy has persisted when comparing model results with trimlines assumed to indicate past maximum ice surface elevations10–16. Numerical models have traditionally overestimated the AIF’s thickness by between 400–1000 metres during the local LGM20,22, thus misrepresenting the supposedly thin and topographically constrained nature of the ice field11,18. This ice thickness overestimation has been attributed to oversimplified ice flow dynamics and the coarse spatial resolutions of previous AIF-wide models, computationally restricted to 1–2 km22. However, the discrepancy also caused some studies20,28,29 to question whether trimlines in the Alps are true indicators of maximum ice surface elevations, or rather of warm-to-cold-based ice transitions, as was previously debated in other glaciated regions, such as Britain and Scandinavia for instance30.

To address this conundrum, we apply a different approach by modelling the AIF with the Instructed Glacier Model (IGM)31,32, a thermo-mechanically coupled three-dimensional glacier evolution model. IGM makes use of recent improvements in physics-informed machine learning to accelerate the high-order Blatter-Pattyn solver of ice flow33, enabling efficient computation on graphics processing units (GPU). The resulting model simulates glacier motion with fidelity while offering the same ability to perturbate parameters controlling ice properties as traditional central-processing-unit-based models, but at a fraction of the computational cost. Here, we use IGM to simulate the AIF’s transient evolution from 35 to 18 ka, thus bracketing the full LGM period, at a spatial resolution of 300 m, an order of magnitude higher than previously achieved20. Computationally unfeasible with traditional glacier models, our perturbed parameter ensemble experiment shows that increasing spatial resolution substantially reduces model-data misfits in both LGM ice extent and thickness across the Alps (Fig. 1), resulting in a thinner ice field than previous 1–2 km resolution models20,22. This approach enables us to produce a high-resolution, data-consistent, yet physics-based transient simulation of the AIF’s last glaciation. The results have wider implications for addressing open research questions on Quaternary glacial erosion processes, post-glacial isostatic rebound, and paleoclimate in the European Alps. Our study also reveals that GPU-based glacier evolution modelling is a promising tool for better reconstructing other paleo ice fields, but also for past and future ice-sheet-scale modelling. Indeed, by substantially reducing computational costs, this approach permits advances in modelling resolution and parameter space exploration, both essential to resolve complex topographies and model more accurate past and future glaciated environments.

Results
High resolution model setup
To model the transient evolution of the AIF during the LGM, we design a high-resolution (300 m) model setup with IGM31,32 that builds on the Parallel Ice Sheet Model (PISM25) experimental setup of Jouvet et al.22. We implement an equivalent enthalpy module formulated by Aschwanden et al.34 for modelling polythermal ice (see ‘Methods’ section) and use the same input climate fields27 and geothermal heat flux data35. Surface mass balance is computed using a comparable positive degree-days scheme based on Calov & Greve36. We use a pseudo-plastic power slip law37 and parameterize the space-dependent yield stress to be controlled by the effective pressure and frictional strength of basal till38 (see ‘Methods’ section). Glacial isostasy is modelled by coupling IGM with the gFlex lithospheric flexure model39, after Mey et al.21.

The experiment setup of Jouvet et al.22 is however improved by using an Alps-specific climate record (instead of the EPICA40 ice core record) as input to our glacial index scheme22,41 (Fig. 2e). The applied signal 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). Further improvements of the setup include an avalanche scheme re-distributing snow and ice accumulated on steep slopes (>45°) down to the glacier surface, and an elevation-dependent parameterization of the basal till friction angle, enabling a spatially-variable bed strength (see ‘Methods’ section). In the main overdeepened valleys of the Alps (e.g. the Rhône valley), the thickness of unconsolidated infill sediments present before the last glaciation remains debated21,44. Here, we test the impact of valley-fill sediments on the ice field’s LGM geometry by using two distinct digital elevation models within our ensemble (see Supplementary Table 1, ‘Methods’ section). In the first digital elevation model (from Jouvet et al.22), only present-day lakes and glaciers are removed45, while the second (from Mey et al.21) additionally removes all valley-fill sediments from major valleys across the Alps.

In this study, IGM makes use of a physics-informed convolutional neural network to emulate a high-order (three-dimensional) Blatter-Pattyn ice-flow solver32 (see ‘Methods’ section). This methodology reduces computational costs by several orders of magnitude through efficient parallelization on GPU. This approach enables us to run AIF-wide simulations (542,700 km2) at a spatial resolution of 300 m (grid size: 3006x2006). Over the entire Alps, a 300 m grid increases the total cell number for which to compute ice flow by a factor of ~11 relative to 1 km (best achieved so far20). While running a 20 kyr-long simulation at 300 m would require ~2.5 years with PISM using a 32-core (3.7 GHz) computer (or ~6 months using a state-of-the-art high-performance computer), it is achieved in only ~2.5 days with IGM using a single GPU (Nvidia RTX 4090).

Taking advantage of such capabilities, and after validating IGM’s suitability through quantitative comparison with earlier PISM results (see ‘Methods’ section), we ran AIF-wide simulations covering the full LGM period, from 35 to 18 ka. To assess model-data agreement (see below) and investigate the model sensitivity to the choice of parameters for a wider range of ice field evolution scenarios, we perform 100 simulations. Within this ensemble, we vary 10 key model parameters drawn from eight different components of the glacier system which modify both ice properties and boundary conditions (see Supplementary Table 1). The 100 parameter sets are sampled between given parameter ranges using a Latin hypercube algorithm46, optimised with a maximin criterion. This commonly used technique provides efficient and homogeneous sampling of parameter spaces with high dimensionality, helping to optimally explore the variety of possible model responses46.

Model-data comparison scheme
Model-data agreement is first evaluated for ice extent by comparing modelled maximum ice thickness fields against an empirical reconstruction of AIF margins at the LGM originally produced by Ehlers et al.17. In this study, we use the newest version of this outline (Fig. 1a), which was updated by a series of studies5,8,47–53 providing more constraints on LGM margins of main AIF outlet glaciers (see ‘Methods’ section), and which is used widely within the community5,6,8,18,52,54. The fit for each ensemble simulation is quantitatively assessed by computing to total number of overlapping pixels between the modelled LGM ice margin and rasterized polygons mapped around the empirical LGM outline (see Figs. 1a, 2b, Supplementary Fig. 19). For this test, LGM moraine polygons are mapped exclusively in regions of high confidence in the position and broad LGM timing (35–18 kyr BP) of the AIF margin due to abundant mapping and dating5,6,18. The overlapping pixel counts are subsequently normalized to obtain a score between 0 (worse-fit simulation) and 1 (best-fit simulation). As a first model-data comparisons sieve, any simulations scoring below 0.8 (n=58) at this test were excluded (Fig. 2b). Here, the modelled LGM state for each simulation is obtained by computing the maximum thickness reached in each pixel at any time during the simulation, thus providing a time-independent map of maximum modelled thickness and extent (Fig. 1a). The precise dated timing of the LGM in specific regions is thus not directly compared with the timing of the modelled LGM as part of this model-data comparisons scheme, which exclusively aims at finding simulations with the most data-consistent AIF geometry (see Supplementary Fig. 19).

Secondly, we assess model-data fit in ice thickness for the remaining simulations by computing the difference between modelled maximum ice surface elevations and 353 reported trimline elevations from the literature10–16 (Figs. 1b, 2d, see ‘Methods’ section). To ensure a direct comparison, modelled surface elevations are corrected by computing and adding the LGM lithospheric deflection at each trimline location. Finding simulations with lowest mean and standard deviation values in elevation difference helped minimise both the overall model-data misfit in ice thickness and its scatter, the latter acting as a proxy for spatial variability in model-data agreement (Fig. 2d, f). Simulations with mean and standard deviation values exceeding 160 m and 100 m elevation difference, respectively, were sieved out (n=32). This quantitative three-sieve model-data comparison procedure was designed to isolate a group of 10 best-fit simulations (10% of the ensemble) for further inspection. Finally, a visual check of these 10 simulations by three experts enabled to identify two outlier simulations which produce far too extensive ice margins throughout the Western and Northwestern Alps. After removing these two outliers, the final eight Not-Ruled-Out-Yet simulations (NROYs) are found to be indistinguishable in quality of ice-extent fit with the empirical LGM outline and are thus used for all subsequent quantitative results reported in this study, using NROY means and standard deviations in given metrics (Fig. 2). For visualisation and discussion purposes only, simulation 37 was further selected as our best-fit simulation, due to presenting the best combination of extent and thickness model-data fit.

Improved LGM model-data fit
Overall, our model produces a thinner LGM ice field resulting in more valley-confined ice flow (Fig. 3), and leading to a substantial increase in model-data fit in ice surface elevation. Using the same trimline dataset10 as Jouvet et al.22 and Seguinot et al.20 (n=175), the mean model-data misfit with observed trimlines for our best-fit IGM simulations is +191±15 m (NROY mean). The positive bias in modelled ice thickness is reduced by 200% and 450% compared with these earlier studies, respectively. When using a more extensive (n=353), spatially widespread, and quality-controlled trimline dataset10,11,13–16 (Figs. 1, 2, 4b, see ‘Methods’ section), the mean misfit value is further reduced to +146±12 m (NROY mean). Across the Alps, modelled maximum ice surface elevations are positively and well correlated with reported trimline elevations (NROY mean r2=0.92) (see Supplementary Figs. 2, 3), with misfit values for all trimlines presenting limited scatter (NROY-mean stdev=94±4 m). Despite varying 10 key parameters conservatively (see Supplementary Table 1), all 100 simulations present mean trimline misfit values that remain below previous and coarser-resolution models’ best estimates20,22. This implies that increasing model resolution from 2 km to 300 m systematically produces a thinner LGM AIF throughout the Alps (Fig. 3a), on average by 218±15 m (NROY mean) relative to Jouvet et al.22. Indeed, maximum AIF volume (57.6±2.1x103 km3) is here ~35% lower than when modelled at 2 km (88.3x103 km3)22, while maximum AIF area is <2% different (Fig. 3). The resulting ice thickness difference is spatially heterogeneous, however, and most prominent for the Rhein and Garda outlet glacier catchments. The NROYs, for instance, model both the adjacent Rhein and Linth outlet glaciers reaching the mapped LGM moraines (Figs. 1, 3) while keeping the Hörnli nunatak (47.4°N, 8.95°E) and its surrounding hills ice free. Although well documented7,55, this LGM ice-field configuration was never achieved by previous AIF-wide simulations, which instead modelled ice covering the entire region20,22.

The large and recurring model-data thickness offset obtained in previous modelling efforts caused certain studies20,28,29 to suggest that trimlines in the Alps may represent a warm-to-cold-based ice transition. Across our 100 simulations, we find highly variable modelled basal ice temperatures at the locations of trimlines during the LGM, with no apparent clusters towards distinct temperatures. This suggests no relationship exists in our modelling between thermal boundaries of basal ice and trimline formation (see Supplementary Fig. 4). These results suggest that in the Alps, observed trimlines are likely true indicators of maximum ice surface elevations during Quaternary glacial maxima.

Motivated by potential uncertainties in the input climate forcing27 and our dependence to the glacial index approach41, our ensemble features a catchment-specific precipitation offset scheme (see ‘Methods’ section). The latter permits simulation-specific climate variations with possible localised changes in input precipitation of between −40% and +50%. Together with higher-resolution modelling, our approach also improves model-data fit in LGM ice extent throughout the AIF, relative to previous work (Figs. 1a, 4a). Where Jouvet et al.22’s PISM simulation produced AIF margins more extensive than empirical reconstructions (i.e. in the Rhein, Reuss, Durance, Jura, Vittorio-Veneto, and Astico regions), our best-fit simulations produce less extensive, more data-consistent ice margins. Conversely, the fit is improved where previously-modelled AIF margins22 did not reach the mapped ice margins (i.e. in the Drau, Inn, Isar, and Lyon regions) (Fig. 4a). The LGM extent of the Drau outlet glacier (Figs. 1a, 4a), for instance, was notoriously challenging to model accurately20,22,24. Our NROYs reproduce its extent within 10 km of well-preserved and mapped local terminal moraines9, and present a margin shape consistent with the consensus empirical outline17. However, we note three instances, i.e. the Rhône, Dora-Baltea and Isère outlet glaciers, where our NROYs model too extensive ice margins leading to a worse fit with data than Jouvet et al.22’s model (Fig. 4a). We believe remaining model-data misfits in ice thickness (+146±12 m, NROY mean) and in the extent of certain outlet glaciers can either be, 1) Reduced further with higher-than-300 m resolution simulations and more extensive parameter-space explorations using larger ensembles, 2) Related to unavoidable uncertainties inherent to paleotrimline and LGM moraine identification due to challenges with landform preservation and dating, and 3) Related to uncertainties and biases in our input climate and surface mass balance parameterisation. A combination of impact from these three mechanisms is also considered likely. For the latter, future improvements in the input climate through higher-resolution regional climate modelling and more transient simulations accompanied by a more complex energy balance model may help reduce some of the remaining model-data misfits.

Intra-NROY variability
Overall, our NROYs show few clusters within the sampled ranges of the 10 ensemble-varying parameters (see Supplementary Fig. 5). An exception is the Glen’s flow law enhancement factor (E), with NROY values all greater than one. Thus, good model-data fit may require more deformable ice than modelled with default flow law constants (E=1, n=3). The lack of clusters in other parameters indicates a better model-data fit can be obtained with highly variable parameter values, making optimal parameter configurations difficult to predict. This variability justifies a perturbed ensemble approach to maximise model-data fit, as numerous local minima may exist. Moreover, two out of eight NROYs were achieved using a basal topography without valley-fill sediments, while the remaining six were not (see Supplementary Fig. 5). A sensitivity analysis on the best-fit simulation (37) further reveals that at the AIF scale, the removal of thick (>300 m) valley-fill sediments prior to glacier advance does not significantly impact ice flux, glacier build-up rates, ice-field geometry, and thus model-data fit at the LGM (see Supplementary Fig. 6). Therefore, we find inaccurate thickness estimations of valley-fill deposits prior to the last glaciation have little impact on the quality of AIF-wide models of the LGM.

Methods
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:

∂h/∂t + ∇·(u·h) = MB, (1)

where ∇ represents the divergence operator for the horizontal variables (x,y), u· = (u·,v·) denotes the vertically-averaged horizontal ice velocity field, and 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:

D = Aτn (2)

where D and τ 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 τb to the sliding velocity ub as follows:

τb = -c|ub|m-1ub, (3)

where m>0, c=c(x,y)>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 formulation of ice flow and temperature, as shown in Aschwanden et al.34. Ice enthalpy, E, is a function of temperature, T, and water content, ω:

E(T,ω,p) = cp(T + Lω/cp), if T < Tmp(p), ω = 0, (4)

where the temperature, Tmp, and enthalpy, Emp, at pressure-melting point of ice are defined by

Tmp = T0 - βp, (5)

Emp = ci Tmp(p) - Tref (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(T,ω) = Aexp(-Q/(RTpa)) (1+181.25ω), (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, Wtill, is computed from the basal melt rate, mb, obtained from the enthalpy as follows:

∂Wtill/∂z = mb/ρw - Cdr, (8)

where Cdr is a simple drainage parameter. The till layer is assumed to be saturated when the basal water layer thickness reaches a cap value of Wtillmax = 2 m. The effective thickness of water within the till layer Ntill is computed from the saturation ratio s=Wtill/Wtillmax by the formula78:

Ntill = min(p, N0 (1-s)/(1-s+C0)), (9)

where C0 is a constant. The ice solver in (3) is defined by the parameters of the Coulomb law56 that involves the effective pressure in the till Ntill:

τb = Ntill tan(ϕ) ub/|ub|, (10)

where ϕ is the till friction angle25. Here, using the assumption that basal materials are generally weaker (softer sediments) in valley troughs than over mountaintops38, ϕ is parameterised to be a piece-wise linear function of bed elevation, b:

ϕ(b) = ϕmin, b ≤ bmin; ϕmin + (b-bmin) (ϕmax-ϕmin)/(bmax-bmin), bmin < b < bmax; ϕmax, b ≥ bmax (11)

Between ensemble simulations, we modify the elevation-dependency of ϕ by varying upper and lower bed elevation thresholds (bmin, bmax) between ranges of (−500, −100) metres, and (2400, 3000) metres, respectively, while till friction angle thresholds (ϕmin, ϕ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 ubth 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 work-flow 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 the spatially-variable 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) reveal 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–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 outlines 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.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of Leger et al: a data consistent model of the last glaciation in the Alps achieved with physics-driven AI\nThank you for the opportunity to review this paper. This manuscript seeks to improve the modelling the European Icefield\nusing the Instructed Glacier Model. The use of GPUs means that processing is fast and efficient, allowing higher spatial\nresolution across the domain. The ability to model icefield-wide at 300m resolution is an extraordinary achievement. This is\nhigh enough to capture the deep, narrow alpine valleys, a major control on ice flow and ice dynamics that are typically\nneglected by more computer intensive models, and the use of higher order Blatter-Pattyn physics ensures that the most\nimportant stresses are accounted for. The model-data comparison framework to score the ensemble is state of the art. The\nresults of the paper and especially the methodological steps forward gained have far-reaching implications.\nThe surface mass balance scheme used is pretty simple however and perhaps this is the next area that the authors could\ntackle to improve!\nThis manuscript is an exciting paper that I enjoyed reading. It is truly transformative science at its best and a major leap\nforwards. My congratulations to the authors. I have only a few suggestions, mostly for clarity, and to help the paper fulfil its\npotential as a highly cited transformational piece of work.\nMostly these comments are to do with: design of the model-data comparison; and the wider implications in the discussion.\nModel data comparison:\nThe empirical outline used in figure 1 is used somewhat uncritically and I think more discussion is needed (perhaps in the\nMethods or Supplement) more analysis of this limit. What is the uncertainty in the limit? What ages are the moraines and are\nthey all the same age? When you score your model based on the best fit to moraines, this assumes that all moraines are the\nsame age and are well dated. However, we know that topography has a significant influence on moraine formation. Where\nyou have model underfit (e.g. GA glacier), what is the confidence in the age of the moraines?\nThe Ehlers limit is fairly dated now and on reading it wasn’t clear to me if this was used as published or if it was updated\nbased on the more recent works cited. This needs to be clearer as although trimlines are used largely to score the ensemble\nthe moraines are used to find the best fit scenario. How well dated are the moraines and the trimlines?\nTrimlines are also challenging to date – often boulders and erratics here will be transported supraglacially and therefore\nhave an inheritance when used in cosmogenic exposure age dating. In the UK for example, trimlines have proven hard to\ndate and age is debated.\nOverall, some more discussion somewhere about the data used, the ages of the moraines, and how much confidence there\nis in it would be helpful. This could be in the Methods and SI. I also wonder if this may help explain where you have\nunder/overfit.\nThe ‘quality controlled’ reference to trimlines in caption of figure1 – this was not clear to me what this entailed.\nWider implications:\nI wondered about the wider implications of the work and think that the authors could do more agenda setting in the\ndiscussion. The Discussion is rather inward looking and the wider implications to a wider audience can come across more\nclearly.\nThe use of trimlines to realistically constrain the model ensemble members is exciting and I can see this spawning a new\ngeneration of mountaineer glaciologists climbing peaks to obtain cosmogenic nuclide ages. If the use of trimlines rather than\nmoraines (more commonly used) to evaluate the ensemble makes for a systematically better product, then you can highlight\nthis to glacial geologists and indicate the kinds of data needed for model-data comparisons. I think that this would bring the\nwork to the forefront of a wider audience and would increase readership and applicability beyond the modelling community.\nEmphasise the need for more well-dated, well constrained trimline data in those parts of the world with currently rapidly\nchanging icefields (Canada, Alaska, Patagonia, Himalaya, Antarctic peninsula), where the impacts of glacier recession will\nbe significant, and highlight that constraining past rates of thinning would be especially beneficial to improving future\nprojections. The kinds of data needed could be more clearly signposted.\nThe perturbed ensemble approach to maximise model-data fit will be critical in forward modelling these icefields and\ndemonstrates the importance of hindcasts and data model integration to calibrate projections. This can be emphasised more\n– there is a rather vague reference to Greenland and antarctica in the conclusion but here, the thicker icesheets have fewer\nnunataks and it is likely the smaller mountain icefields more similar to the Alps that this approach will be most beneficial.\nIs the model set up useful for forward modelling of the Alps ?\nMinor comments\nI wondered how lakes that formed in over deepenings were accounted for and if frontal ablation was simply neglected. Does\nthe model assume that all glaciers are land terminating? Frontal ablation is likely to be a significant driver of ice dynamics so\nalthough this is not an issue to be corrected it would warrant a minor sentence or two to clarify.\nIn the experiment design (line 218) please clarify: are the parameter spaces weighted or selected according to probability or\nknown ranges?\nLine 4 -7: It is helpful to cite some references for mentioned ‘an abundant library of mapped and dated ice-contact glacioenic\nlandforms and sediments separately for reader’s interest.\nLine 15-17: Name the glacier evolution models (for example with some references)\nLine 39: What is the transient time-step of model simulation results, is it the same as the temporal resolution of climate\ndataset?\nLine 169: (I would expand more on ensemble for clarification: is it ensemble of 100 parameterized model simulations\nresulting from a group of parameter sets and /or a group of climate input used in the study)\nhow is the initial range for each parameter selected? Is it from the output of the model simulations or based on the model\ndefault values or from literatures?\nLine 225 -227: Any explanation for why the maximum AIF areas outside the range of 156-165X 103 produce substantial\nunder-or overestimations of the LGM ice extent reconstruction?\nLine 228- 230: Could it be because of the deep learning algorithm that the model is producing similar patterns for outlet\nglaciers despite extensive changes in the parameters? does this also mean these outlet glaciers control the evolution and\nbehaviour of the entire AIF?\nLine 241-243: how is this natural break for mean and standard deviation i.e. 180 m and 100 m in elevation difference\nchosen? (is it from sieved simulation data set?)\nLine 247: the best fit to moraines ensemble member selection: was how was this quantified? Visually?\nLine 280 – unprecedented is a very strong word, can it be justified? I would avoid using more superfluous language\npersonally.\nLine 326-327: why is that so? Is it because this outlet glacier is influenced by more local phenomena than the Alps-wide\nphenomenon and the ensemble parameter sets could not capture the processes and drivers?\nLine 408: explain briefly what high-magnitude venturi effects mean (not necessarily all glaciologists know what it means)\nLine 557 – how good is the chronological evidence here?\nLine 637 – how else could it inform future investigations – see comments above.\nLione 647 – there are also implications for data as well\nLine 656 – this is very vague and I suggest rewriting this paragraph. Greenland & WAIS/EAIS have very few nunataks and\nrequires a different approach. Can a very small number of spatially limited dispersed nunataks improve modelling here?\nLine 659 – a rather superficial sentence to end on. You can finish with a more exciting call to arms and set the agenda for the\ncommunity here.\nSome limitations and how they can be improvised have been explained in the discussion sections (Line 630- 637).\nHowever, it would be beneficial to clarify the major limitations of such high-resolution, 3D -AI-based models and outline the\nkey considerations when integrating multiple physical processes into a single, high-functioning model.\nDiscussion\nThis section is well-explained and effectively discussed in relation to previous research. It offers valuable new insights into\nboth glacial dynamics and landscape evolution. The revelation of localized ice flow speedups in narrow valleys and\ntopographic bottlenecks, along with rapid valley deepening in areas like the Rhône valley, highlights the critical role of high-\nresolution, physics-informed models in addressing unresolved questions in large-scale past glacier modelling. Additionally,\nthe inclusion of future research perspectives and avenues with the IGM model is appreciated. This will be more useful for the\ncryosphere community and quaternary researchers.\nMethods\nMore information on model ensemble evaluation would be useful here.\nMore information on the data used for model evaluation would also be useful.\nCan IGM be advised to use for reconstructing a single glacier elsewhere in the world with same ensembled parameter sets?)\nBethan Davies", "reviewer_2": "Review of Leger et al\nThis study presents a major breakthrough in glacier modeling capability. New information and understanding gleaned from\npaleo-glacier simulations, and future glacier simulations, is most robust when models run in ensemble mode, that is, large\nnumbers of simulations that allow one to account for a range in model inputs. Ensemble results also provide a way to\ncharacterize the uncertainty of numerical simulations. But running ensembles does not occur frequently because most\nnumerical glacier models are too computationally demanding. This study relies on emulators and other marine learning\ntechniques to speed up the modeling process, and allows for the ensemble approach.\nThe authors test this model capability in an oft-used target for glacier modelling (the European Alps) and builds on recent\nwork that provides context for the new cutting edge approach. This nicely described and illustrated study has me convinced\nthat it is deserving of publication with minor revisions.\nI have some comments more significant than others, these I keep in the line-by-line section below, a summary here:\n- Better explain choice in sieve ranges\n- Use 8 NROYs to quantify uncertainty of final result\nLine by line\nAbstract (no line numbers) just some word-smithing suggestions, results ‘have’ implications instead of ‘yield’ implications.\nSuggest saying something about how the decreased computation time allows users to model glaciers with ensembles, the\nway to go to test parametrizations, allows for uncertainty quantification, the list goes on. This goes well beyond modeling\nglaciers in the Alps better.\nLine 7, ‘one of the better reconstructed’ reword..\nLine 12, ‘also help reconstruct’ reword…\nLine 28 – these trimlines are debated all over the world (not just Britain and Scandinavia, but of course you can’t cite all)\nLine 31 – this paragraph sure has a lot of acronyms. Acronyms are overrated (!), do we really need them at all anymore? Are\nthere space constraints and word counts? Suggest doing away with them, it just makes writing more inaccessible and\ndifficult to understand.\nLine 47, this reference of 1-2 km grid should come several sentences earlier when discussing that older model\nLine 58, upon reflection, half or more than half of the introduction is just a summary of the results and the paper – that’s not\nreally introduction. Yes, many times introductions close with perhaps a final paragraph on ‘what we did’, but this seems a bit\nmuch. Suggest reducing summary text that appears in the introduction, and in place adding more actual intro material. As is,\nit’s like reading results, then going to the results section to read again the same results..\nLine 107, not sure I understand the “Are shown; blah” syntax… Change?\nLine 143, are you sure “results” is the best description of these first paragraphs – seems more like methods to me. Can you\ncome up with a title that’s neither, like “the Instructed Glacier Model” or something?\nLine 149, ‘follows and improves’ well I hope it would improve! It sort of goes without saying, can you replace with ‘builds on’\n?\nLine 156, need word ‘materials’?\nLine 166, In ‘the’ main overdeepened valleys…\nLine 190 – figure 2/study design… I recommend explaining in the text where the sieve boundaries come from. How are they\nderived? Why not choose a narrower range? Or a broader range? Is it arbitrary, or based on something? Sorry if I missed it.\nLine 203, can you add a comma in that big number (6 digits)\nLine 207, can you do away with the footnote and just add the sentence to the paragraph?\nLine 212, recommend instead of ‘we run a series’ write ‘we ran 100 simulations from 35 to 18 ka’. I think at this point we\nknow it’s using the 300 m grid.\nLine 213, data-model agreement is mentioned, but I, for one, didn’t understand yet what data – were they introduced yet, or\nmentioned yet?\nLine 226 – how is this range chosen?\nLine 233, might be nice here to add that 18 sims are left\nLine 244, say ‘3-step’ but that confused me because text describes it as a 2-step. I realize there is a volume step that is\nderived from ice-thickness (?) so it’s kind of like there are two parts of a second step – but I found the description of this as a\n‘3-step’ process a little confusing. Maybe wording will help.\nLine 246, simulation 37. I agree with your choice to pick one for illustration. I got hung up on this a little bit, because if, in the\nend, you just pick one among the 100 that matches the data best, then why do the 3-step sieve process? Why have 8\nNROYs and not just pick the best? I think you have 8 NROYs (not just one) so you can have some attempt at assigning an\nuncertainty to your result, no? Maybe you can articulate this better, if that’s the case. Otherwise, ignore the sieving and just\nchoose the single best.\nLine 274, ‘thus permitting to compute’ reword..\nLine 283, reword for clarity\nLine 293, suggest phrasing ‘… increasing model resolution from x km to x m systematically …’\nLine 297, resulting ‘thickness difference’ ?\nLine 302, replace ‘modelled’ with ‘achieved’\nLine 308, remove comma\nLine 308/309, huh? Reword?\nLine 329, However, we…\nLine 331, ‘due to too extensive…’ reword\nLine 402, are these two phenomena not just one?\nLine 470, suggest ‘cold based regions cover 14±3% of the AIF (NROY mean)(Fig. 6b), more than the 11% modelled at 2 km\nresolution’\nLines 695 and 890, find a better word/phrase to replace ‘hardly’ with, which is a bit colloquial.\nJason Briner", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_4": "In `A data-consistent model of the last glaciation in the Alps achieved with physics-driven AI' the authors present a much\nimproved simulation of the Pleistocene glaciation of the Alps, which is facilitated through the use of an efficient ice physics\nemulator. Leveraging the fact that the emulator allows for fast simulations, the authors produce an ensemble of randomly\nvaried simulations. 8 of these produce ice geometries that are in substantially better agreement with observations of ice\nextent and especially maximum elevation (as deduced from trimlines) than previous studies. Based on the resulting models,\nthe authors draw a number of interesting conclusions about the climatic and glaciological conditions that must have formed\nthe Alpine ice field.\nI find this manuscript to be both very interesting and of exceptional clarity. I have no major issues with the work, which I find\nto be methodologically sound as well as a demonstration of an exciting new hybrid of machine learning and classic\nglaciological approaches. Indeed, please do not take my relatively few detailed comments (which are below) to indicate that\nI did not read the manuscript carefully - I read it several times and found that previous points or questions that I had noted\nwere already answered clearly elsewhere in the text. I also believe that this manuscript will be of significant interest to a\nvariety of audiences.\nL34: \"so-called high-order (three-dimensional)\". I'm not sure that any of these terms are relevant to a general audience -\nthese terms (or the subtexts associated with them) are terms-of-art for technical glacier modelers and may not need to\nappear here.\nL37: What does \"identical levels of parameter control\" mean?\nL159: Awkward sentence.\nL146--247: I appreciate that Nature-type journals don't want methods in the manuscript, which is why these sections (which\nclearly are methods) are labeled as `results'. But maybe it's okay to just call them methods anyways? Or have an\n`experimental setup' section before results?\nL241: Slightly more detail on what constitutes a `natural break' and how robust that selection is would be helpful for\nreproducibility.\nL279: What does `topographically constrained' mean?\nL282: It seems that the model remains considerably biased relative to observations, improvement versus previous works\nnotwithstanding. I think it would be of value to try to explain why this remaining bias exists - especially relative to mean ice\nthickness or something - and if it's something that needs to be considered when contextualizing results.\nL408: Why would it be the case that the inclusion of longitudinal stresses produces higher velocities?\nL676: What is the numerical method used to update the ice geometry? My understanding is that IGM only emulates the\nvelocity field. If that's not the case, then how is it ensured that essential conservation laws (especially mass) are upheld?"}, "author_response": "Modifications related to main comments by reviewer 1\n######################################################################\nModel-data comparison\n- We thank the reviewer for raising this very good and important point here on the model-\ndata comparison scheme. Please find below some more information about the original\napproach, but also regarding a new, more quantitative approach.\n- The original model data comparison approach we followed included an initial step\n(sieve 1) that was semi-quantitative, as based on the visual observation of AIF-wide\nmodel-data fit in ice extent with the empirical LGM outline. This was done for every\nensemble simulation, as the experiment ran, and was noted down in a spreadsheet. This\nled to the observation that any maximum modelled AIF areal extent above and below\nspecific capping numbers resulted in clearly too large or too retreated ice extent,\nand this homogeneously through the AIF (see reply to other reviewer 1 comment regarding\nwhy we think that is). The next two sieves of model-data comparison, however, were fully\nquantitative, using the trimline data with the mean offset, and stdev of offset.\nThis approach resulted in 8 NROYS that appear to be the best-fitting ensemble\nsimulations, and that are indistinguishable in model-data fit quality in terms of ice\nextent (when viewing them on a map) whilst minimizing the misfit in ice surface\nelevation with trimlines. However, we agree that a more quantitative first step would be\npreferable: especially since it was responsible for a large number of simulations being\ndiscarded (82 % originally).\nThus, in response to the comment/question, we designed and tested a different, more\nquantitative model-data comparison and sieving workflow: presented here:\n1) We fist rasterize the empirical LGM outline shapefile from Ehlers (and many other\nstudies, see reply to reviewer 1 comment on this) using a polygon buffer wide enough to\nbe rendered in a 300 m resolution grid. Here, after conducting a literature review, we only\nkeep the margins that are described as high confidence due to good moraine\npreservation and robust dating attributing the margin to around the LGM (anytime\nbetween 35 and 18 kyr BP works for us as our model-data comparison is time-\nindependent within that broad LGM timeframe). We do not use the entire LGM outline:\nwhich yields high uncertainty in certain regions (due to lack of precise mapping and\ndating). We make sure to have data around the full periphery in a broadly homogeneous\nway, though. This literature review was conducted as part of the ALPICE database\nproject (in prep) by Dr. Sarah Kamleitner.\n2) For each ensemble simulation, we isolate the modelled LGM ice margin by using\ntime-independent maximum ice thickness reached for each pixel throughout the full\nsimulation time (e.g. figure 1A). We isolate the ice margin using raster dilation between\nmodelled ice and no-ice masks (using a kernel of 7x7 pixels). This enables to\nautomatically isolate the modelled LGM ice margin as a rasterized buffer zone (see new\nsupplementary Figure 19 and caption).\n3) We then add the rasterized empirical LGM outline and the modelled LGM ice\nmargins. We then calculate the total number of pixels that overlap between the two. The\nhigher that number, the better the fit in ice extent for the LGM. We calculate that number\nfor all ensemble simulations, normalize the full list between 0 (worse sim) and 1 (best\nsim). This gives us a final score between 0 and 1 for each ensemble simulation. This\nunique number estimates the fit in LGM ice extent across the entire AIF.\nIf we apply this method and use this score as our first sieve: by only keeping simulations\nwith a score above 0.8, for instance, and then follow with the two other sieves\ncomparing modelled max surface elevation with trimline data, we find that we end up\nisolating the same NROYs than our 8 original NROYs, using a fully quantitative approach.\nHowever, this is not perfect and we find this also produces two extra outlier simulations\nbesides the 8 NROYS, which are not filtered by this sieving procedure, and which clearly\nrepresent a worse fit in ice extent than the 8 indistinguishable NROYs, when their map of\nmaximum time-independent ice thickness is visualized (e.g. Fig. 1a). In these 2 outliers,\nwe find modelled maximum ice extent in the Western and Northwestern AIF (which we\nknow to be well dated and mapped) ends up being far too extensive. We believe this is\ndue to the human eye using many criteria for an assessment of model data fit in ice\nextent (e.g. overall fit and magnitude of overlap, underlap, with implicit weights), while\nsuch a complex assessment is very difficult to fully automatize, and in our fully\nquantitative method, only one hard criteria is used. Increasing the harshness of the\nsieves further ends up getting rid of some of the 8 NROYs, whilst keeping 1 or 2 of the\noutliers. If kept like this, the resulting \"new\" NROYS (n =10) would become\ndistinguishable in quality: one could clearly see that some simulations within the\nNROYS are better than others in model-data fit. This would be a bit problematic.\nTherefore, we propose this new approach (preferred):\n1 - Our new quantitative scoring of ice extent fit (presented above) is used: and a sieve\nof 0.8 is used (simulations scoring below that are excluded). This is sieve number 1.\n2- We fully get rid of the upper and lower cap on maximum total AIF modelled areal\nextent from the procedure: which was semi-quantitative.\n3 - As per the original scheme, the mean and stdev of the misfit on maximum ice\nsurface elevation are used as sieves number 2 and 3. Here, threshold values of 160 m\nand 100 m are used, respectively, in order to isolate 10 best-fit simulations.\n4 - However, we keep a final visual assessment of the 10 remaining simulations to\nassess whether they are truly indistinguishable in model-data fit quality with regards to\nice extent fit the empirical outline.\nDue to the small number of simulations left (n=10), this visual check is relatively quick\nand easy, and is here conducted by 3 different experts within our co-author team.\n- From this final check, 2 clear outliers are identified (in consensus) for their worse fit,\ndue to far too extensive ice being modelled along the Western and Northwestern AIF\nmargins, and are thus removed (Fig. 2f).\n-We prefer this approach as the human visual observation and expert judgement only\ncomes in the end, and filters out only 2 simulations (2 % of the full ensemble). In the\noriginal approach, a semi-qualitative part of the procedure (sieve 1) was also based on\nvisual observations and filtered out 82 % of the ensemble. We thus consider this new\nmethod to be more robust, less subjective, and more reproducible. In line with the\nimplementation of this new procedure, we updated the method and results sections,\nand figure 2 accordingly, and added a supplementary figure with a caption that explains\nthe workflow for scoring model-data fit in ice extent using a quantitative scoring\nscheme.\nTo summarize:\nThanks to the reviewer's comment, we modified the workflow and manuscript to apply a\nnew, more quantitative, model-data comparison method for the LGM ice extent analysis.\nThe results and list of 8 NROYs are unchanged. However, in order to achieve the same\nquality in results as before, we find we still have to apply one visual assessment of the\nfinal pool of simulations that passed all quantitative sieves (n=10). This final check\nresults in the removal of 2 outlier simulations. We consider this new approach\nrepresents a much less impactful human intervention than before (2% instead of 82%)\nand is therefore much less subjective. However, if the editor and/or reviewer prefers the\noriginal approach, we are of course happy to reverse back to that as well.\n#########\nFurthermore, in response to this comment: “The Ehlers limit is fairly dated now and on\nreading it wasn’t clear to me if this was used as published or if it was updated based on\nthe more recent works cited.”\nWe apologise for not making this clear enough, we have now added a sentence after\nintroducing the outline to better explain that the outline used in this study has indeed\nbeen updated by numerous studies since Ehlers et al., and thus represents the most up\nto date level of knowledge (at the AIF scale) on the position of the AIF margins during the\nLGM:\n“Model-data agreement is first evaluated for ice extent by comparing modelled\nmaximum ice thickness fields against an empirical reconstruction of AIF margins at the\nLGM originally produced by Ehlers et al.17. In this study, we use the newest version of this\noutline (Fig. 1a), which was updated by a series of studies5,8,47,49–54 providing new\nconstraints on LGM margins of main AIF outlet glaciers (see ‘Methods’ section), and\nwhich is used widely within the community5,6,8,18,47,48.”\nWe have also added an entire paragraph to the methods section on this: called “The\nempirical LGM outline of the AIF” which tracks the history of updates made to the\noutline by various studies since. We thank the reviewer as this is a great addition to the\nmethods. However note that adding these references causes us to go over the maximum\nnumber of references allowed (n = 70): which we have already been struggling to respect\nin the original manuscript’s version (we had reached the limit). Would the editor be ok\nwith us going over this threshold? : see the new methods paragraph here:\n“The empirical LGM outline of the AIF\nThe empirical LGM outline of the AIF used in this study was originally produced by an\nAlps-wide compilation of geomorphological and geochronological evidence by Ehlers et\nal.17. Since this key study, the outline has been updated by a series of empirical\ninvestigations improving the quality of LGM margin reconstructions in specific sectors of\nthe Alpine foreland. In this work, we use the most up-to-date version of the outline. Here,\nwe provide a summary of the main studies which have updated the LGM AIF outline\nsince Ehlers et al.17. Gianotti et al.50,51 have made updates to LGM margins of the Ivrea\noutlet glacier in the Dora Baltea region. Braakhekke et al.49 have done so for the Orta\nregion, Kamleitner et al.5 for the Verbano region, Ravazzi et al.53 for the Oglio region,\nMonegato et al.48 for the Garda region, Federici et al.47 for the Gesso region, Ribolini et\nal.8 for the Stura region, and Ivy-Ochs et al.52 for the Dora Riparia region.”\n###########\nFurthermore, in response to this comment: “The ‘quality controlled’ reference to\ntrimlines in caption of figure1 – this was not clear to me what this entailed.”,\nWe apologise for the confusion and have now added new text to the figure 1 caption\nguiding the reader to the paragraph in the methods where we explain our quality-check\nprocedure for the trimline data:\n“Panel b displays the model-data agreement in Alpine Ice Field ice thickness shown for\nour best-fit IGM ensemble simulation (number 37), obtained by computing the\ndifference between 353 quality-controlled trimline elevations from literature10–16 (see\nMethods section: “Empirical trimline elevation data”) and maximum modelled ice\nsurface elevations (time-independent).”\n###########\nFurthermore, in line with reviewer 1’s comments on moraine dating uncertainties: we\nhave added a clearer explanation in our “Model-data comparison scheme” section of\nthe results, regarding the fact that our model-data comparison is time-independent:\n“Here, the modelled LGM state for each simulation is obtained by computing the\nmaximum thickness reached in each pixel at any time during the simulation, thus\nproviding a time-independent map of maximum modelled thickness and extent (Fig. 1a).\nThe precise dated timing of the LGM in specific regions is thus not directly compared\nwith the timing of the modelled LGM as part of this model-data comparison scheme,\nwhich exclusively aims at finding simulations with the most data-consistent AIF\ngeometry (see Supplementary Fig. 19).”\nTherefore, in this work, we do not assume that all moraines are the same age and are\nwell dated. In fact, the moraines described as “LGM” in age by the empirical community,\nand summarized in the outline, could have been deposited anytime between 35 and 18\nkyr BP, the start and end dates of our simulation: since their timing does not impact the\nresults from our simulations and model-data comparison.\nHowever, we fully agree that significant uncertainties lie in the location and confidence\nof the LGM outline due to spatially heterogeneous quality of geomorphological and\ngeochronological constraints and due to data gaps in understudied regions. In the new\nscheme, sieve 1 (which automatically scores simulations based on fit with the extent)\nmakes sure to only uses the outline in regions of higher confidence due to more data,\nand based on a literature review (see above). This information has now been added to\nthe caption of the supplementary figure 19 that explains this scheme for scoring.\n##########\nFurthermore, in response to this comment: “Trimlines are also challenging to date –\noften boulders and erratics here will be transported supraglacially and therefore have an\ninheritance when used in cosmogenic exposure age dating. In the UK for example,\ntrimlines have proven hard to date and age is debated.”\nThe trimlines in the Alps used in this study are not dated. Only bedrock transitions from\nice-moulded to frost-shattered and weathered have been observed and the few\nattempts to date them with cosmogenic dating have resulted in very young bedrock ages\n(Ivy-Ochs, personal communication), due to the very high erosion rates of the high\nalpine environments and its steep topographies. These very high erosion rates have\npushed the community to generally assume that most trimlines in the alps are likely\nLGM in age, as bedrock evidence of maximum thickness limits from pre-LGM Quaternary\nglaciations would have been eroded away, and would not be observable today. We agree\nthis assumption yields uncertainties, and we have now added text in line with this to the\nresults section, also in response to a comment on remaining model bias by reviewer 4\n(see our reply there). However, the fit we obtain between modelled maximum ice surface\nelevations and observed trimlines is now very close (only +146 m of offset, and a R2 of\n0.92), and we believe this good fit shows that either: 1) the assumption that trimlines are\nall LGM in age is robust, or 2) if certain observed trimlines were formed by pre-LGM\nglaciations, then the former maximum ice surface elevations at that time were very\nsimilar to those during the LGM.\n#####################################################################\nWider implications (discussion)\nWe thank the reviewer for raising these very interesting and relevant points.\nWe fully agree that the reduced model-data misfits related to this new ability to model at\nsignificantly higher resolutions over longer timescales, and at much lower cost than\nbefore, thanks to GPU-based models and physics-informed machine learning, should\nalso spark motivation to gather more field data on past glacier and icefield\nreconstructions to bridge data gaps and better calibrate our models. We fully agree that\nthis could be mentioned more in our discussion.\nIndeed, using trimline data in complement to moraine data enables to also obtain a\nmore realistic model in ice thickness, however there are many ways and possible\ndatasets to better constrain a model from field data, some of which are better adapted\nto certain regions due to the topographic settings, the nature of the ice cover, and the\ntype of data available to collect, date, and map.\nIn the Alps, for instance, the high level of topographic constraint on the former icefield\nmakes trimlines a good approach as numerous nunataks existed. This is also the case\nfor many other paleo ice fields (e.g. ice fields in Patagonia, Alaska, New Zealand, Turkey,\nhigh mountain Asia and many others). In other regions, with less high topographies and\nless nunataks, model-data comparison and model hindcasting for paleo\nreconstructions can be improved through other means, such as calibrating models\nagainst ice flow direction using lineations mapping (e.g. Scandinavian ice sheet), using\nerratic transport and provenance databases (e.g. British-Irish Ice Sheet), detailed\nmapping and dating of margin retreat patterns and englacial stratigraphic isochrones\n(e.g. Greenland Ice Sheet), nuclide dipstick models and constraints on past thinning\nrates (e.g. Antarctica) etc... In this work we don’t necessarily want to just focus on a\nwider implications message that mentions trimlines and nunataks only, as these are one\ntype of constraint amongst many others: but more on the fact that high-resolution\nmodelling will enable to reduce model-data misfits and should thus motivate more data\ncollection more generally.\nIn line with these thoughts, we have now modified and added text to the relevant\nparagraph of the discussion:\n“The novel approach presented in this study also has implications for improving the\nreconstruction of other Quaternary ice fields and ice sheet histories. Indeed, achieving\nspatial resolutions required for a data-consistent LGM model of the AIF was previously\\unattainable with traditional central-processing-unit-based models. This study\ndemonstrates that physics-driven machine learning and GPU-based processing makes\nsimulating a continental-scale ice field at high (300 m) resolution with a thermo-\nmechanically coupled, three-dimensional, and high-order model, possible. This is\nmoreover achieved at a fraction of traditional models’ computational costs and with\nidentical parameter modularity. We thus argue that physics-informed, GPU-based\nmodels open the door to a new era of more accurate paleo ice field and ice sheet\nmodelling. In particular, the ability to more accurately model the complex dynamics of\ntopographically-constrained large-scale ice fields, such as the AIF during the LGM,\nthrough high-order modelling using higher spatial resolutions and wider ensemble-type\nexplorations of parameter spaces, will likely reduce model-data misfits in many other\nstudied mountain ranges. Investigations that would likely benefit from this novel\nmodelling approach include the reconstructions of former large-scale ice fields during\nQuaternary glaciations in Patagonia, Alaska, high-mountain Asia, Turkey, Georgia,\nCalifornia, Peru, the Pyrenees, New-Zealand, and Iceland, to name only a few. This\nimprovement in modelling capabilities should also spark new motivation to collect new\nfield data to better constrain the former evolution of ice extent (e.g. with more detailed\nterminal and lateral moraine mapping and dating5), ice thickness and thinning rates (e.g.\nwith new trimlines and cosmogenic nuclide dipstick data69), and ice internal flow\ndirection and dynamics (e.g. with lineation mapping and flow set reconstructions70, new\nerratic transport and provenance databases, borehole data and englacial stratigraphic\nreconstructions71,72) of other former glaciers, ice fields and ice-sheets. This has\nimplications for not only new scientific discoveries in paleoglaciology, but also across\nother disciplines studying processes linked to Quaternary glacial history, including\narchaeology, paleoclimatology, paleoecology, and geomorphology.”\n#####\nFurthermore, in reply to this comment: “The perturbed ensemble approach to maximise\nmodel-data fit will be critical in forward modelling these icefields and demonstrates the\nimportance of hindcasts and data model integration to calibrate projections. This can be\nemphasised more – there is a rather vague reference to Greenland and antarctica in the\nconclusion but here, the thicker icesheets have fewer nunataks and it is likely the\nsmaller mountain icefields more similar to the Alps that this approach will be most\nbeneficial. Is the model set up useful for forward modelling of the Alps ?”\nThe idea of our last discussion paragraph was to focus on the benefit of our new GPU-\nbased, high resolution and ensemble modelling approach to future projections of\ncryospheric change. In this final paragraph, we mainly wanted to focus on the Antarctic\nand Greenland ice sheets for three reasons: 1) The future societal implications, 2)\nbecause modelling their present-day states and future evolution is very computationally\ndemanding at high resolution due to their size, and thus there exist a potent\ncomputational bottleneck there as well (which is not the case with most other ice fields\nand mountain glaciers which can already be modelled at high-resolution with traditional\nhigh-order models e.g. Elmer Ice, due to the short timescales) , and 3) also because\npaleo hindcasting of these large continental scale ice sheets has an impact on future\nprojections since they yield considerable memory and inertia from past climate change\naccumulated over millennial timescales, and which can create multi-century delay in\nresponses to climate change and is thus important for societal implications to better\nconstrain. Our approach thus has implications for the paleo initialisation procedures\nwhen modelling the Antarctic and Greenland Ice sheets. Regarding this last argument\n3): The memory mechanism is not so important for mountain glaciers and smaller ice\nfields which do not have enough ice thickness to accumulate such a memory, and for\nwhich paleo hindcasting in this way is not really an advantage. Inversion schemes\n(which are a different type and much larger type of ensemble) of present-day observed\nvelocity and SMB conditions has been shown to do a good job at representing the state\nof these contemporary smaller glaciers and icefields (e.g. Millan et al. 2022).\nIt is true that GPU-based models such as IGM can also provide great improvements to\ninversion schemes (which are expensive to run) for present-day glaciers and ice fields,\nand also enable global glacier forward modelling in 3D and high order, whilst the latter\nhas so far been restricted to SIA flowline modelling (e.g. with OGGM). However this is a\nvery different type of modelling to what we do in this study, and has already been talked\nabout in detail in the publication of Cook et al. (2023) who used IGM to inverse and\nproject the future of all glaciers in the Alps, and who is currently working on\nimplementing the same technique to all glaciers and icefields in the world. We feel it\nwould be a bit outside the scope of our study to mention this in our discussion.\nHowever, we believe our approach could really help improve the modelling of the future\nevolution of major ice sheets (i.e. Greenland and Antarctica). In line with these thoughts\nand also more minor comments by reviewer 1 regarding the last discussion paragraph,\nwe have tried to modify the text and add new sentences to make this last point less\nvague and clearer:\n“We also believe this new technology may permit more accurate simulations of past and\nfuture Greenland- and Antarctic-ice-sheet change, by facilitating ice-sheet-wide\nsimulations at considerably higher resolutions than previously achieved and at lower\ncomputational costs. This would allow for broader explorations of parameter spaces and\nice-sheet evolution scenarios, with implications for more accurately modelling basal\nconditions and ice-ocean interactions towards grounding lines, crucial for better\nprojecting the future response of continental ice sheets to climate change. Moreover,\nreducing Quaternary model-data misfits, as is achieved here with the AIF, may be\nimportant for improving the paleo initialisation procedures of Antarctic and Greenland\nice sheet models, which greatly impact their future projections. Therefore, we expect\nphysics-informed, GPU-based models to bring advances in projecting future global sea\nlevel rise and environmental change associated with the response of contemporary ice\nsheets to climate change.”\nResponse to reviewers’ line by line comments\n######################################################################\nReviewer 1\n-Reviewer:\n“I wondered how lakes that formed in over deepenings were accounted for and if frontal\nablation was simply neglected. Does the model assume that all glaciers are land\nterminating? Frontal ablation is likely to be a significant driver of ice dynamics so\nalthough this is not an issue to be corrected it would warrant a minor sentence or two to\nclarify.”\nAuthors:\nWe thank the reviewer for raising this point. At the moment IGM (which is still fairly new)\ndoes not have a proglacial lake module (although this is planned as a future point of\nmodel development). In this experiment, the model therefore assumes that all glaciers\nare land terminating. Although not perfect, this assumption is likely to not cause a\nmisrepresentation of the modelled AIF during the LGM (the only event we focus on in this\nstudy), since the majority of over-deepened basins susceptible to host a proglacial lake\nwere eventually filled with ice during maximum AIF expansion. The presence of\nproglacial lakes may slow down the rate of frontal advance during buildup a little\n(although this would need testing), but the magnitude of glacier advance (which is what\nwe exclusively test here) is not likely to be significantly influenced, we believe. This\nassumption is furthermore validated by the quality of fit we obtain between our\nmodelled LGM margins and mapped / dated LGM moraines. However, if one wants to\nmodel the post-LGM deglaciation of the AIF from those margins and test that event\nagainst data, then we believe adding proglacial lakes to the model would be very\nimportant, in order to more accurately estimate the rates of retreat. Here, we\nintentionally stopped our LGM simulations before the major deglaciation of the AIF,\nwhich occurred shortly after 18 kyr BP.\nWe agree this needs to be clarified in the text, and have thus added two sentences to the\nmethods “Surface mass balance” section:\n“Note that no proglacial lake module is implemented in this model setup, meaning that\nall ice is assumed to be land-terminating. We consider this assumption to have little\nimpact on the LGM geometry of the AIF since most overdeepened basins of the Alpine\nforeland were eventually ice-filled during maximum glacier advance.”\n-Reviewer:\n“In the experiment design (line 218) please clarify: are the parameter spaces weighted or\nselected according to probability or known ranges?”\nAuthors:\nMany thanks for pointing this out, indeed the parameter space is initially defined by\ngiven ranges for each of the ensemble-varying parameters (n=10). These ranges, for each\nparameter, are provided in Supplementary Table 1.\nWe have now added this to the relevant sentence:\n“Within this ensemble, we vary 10 key model parameters drawn from eight different\ncomponents of the glacier system (see Supplementary Table 1), sampled between given\nranges using a Latin hypercube algorithm40, optimised with a maximin criterion.”\n-Reviewer:\n“Line 4 -7: It is helpful to cite some references for mentioned ‘an abundant library of\nmapped and dated ice-contact glacioenic landforms and sediments separately for\nreader’s interest.”\nAuthors:\nWe fully agree with the reviewer on this point. Unfortunately, we have already reached\nthe limit number of references allowed (n=70) and thus had to make hard choices as to\nwhich papers to cite. The study covers many subjects from different fields and thus\nrequires numerous key references making us reach the number of 70 references quickly.\nIt will be challenging to add new references (although we agree that would be more\nappropriate). However what we can do is add relevant references that are already being\ncited later in the manuscript within these first sentences to better guide the reader\ntowards these studies. Please see the change here: 18 relevant references are now part\nof these first sentences while originally, only 6 were.\n“The extent, thickness, and flow geometry of the European Alpine Ice Field (AIF) during\nQuaternary glaciations has been studied for more than a century1–4. The result is an\nabundant library of mapped and dated ice-contact glaciogenic landforms and\nsediments (e.g. moraines and tills)5–8, glaciofluvial deposits9, and trimlines10–16, making\nthe extent of the AIF during the Last Glacial Maximum (LGM: ~30-19 ka17,18) one of the\nmost well known in the world11,18.”\n-Reviewer:\n“Line 15-17: Name the glacier evolution models (for example with some references)”\nAuthors:\nWe have added these details to the text. To our knowledge, only PISM and iSOSIA have\nbeen used to model the entire AIF during the LGM or other Quaternary maxima. Other\nstudies within the list of references provided (only two others) have used in-house SIA\nmodels with no names:\n“Under these motivations, over the past 20 years, several studies20–24 have produced\nsimulations of AIF glaciations using glacier evolution models such as the Parallel Ice\nSheet Model (PISM25) and the integrated second-order shallow ice approximation model\n(iSOSIA26), with the main objective to closely fit moraine and trimline evidence.”\n-Reviewer:\n“Line 39: What is the transient time-step of model simulation results, is it the same as\nthe temporal resolution of climate dataset?”\nAuthors:\nThank you for raising this important question. The answer depends on which\ncomponents of the model we look at.\n-The time step of the iceflow module which updates the ice thickness is flexible\ndepending on minimizing computational cost while respecting the CFL condition and\nthus is based on maximum ice velocities modelled in previous time step. In our case,\nduring most of the simulation, it is down to 0.001 yr: ~half a day. This is necessary to\nensure mass conservation and is the case in most glacier evolution models.\n-The time step of the SMB module, i.e. how often the climate fields are updated to\ncompute the annual SMB, using the input climate data (only two snapshots in time: the\nLGM and the PI), the PDD model, and the glacial index scheme: in our case is updated\nmonthly. The annual SMB is then transmitted to the ice-thickness update once a year\n(on 1st of October: end of melt season).\n-The frequency of saved model outputs for post-processing analyses: this is saved every\n50 years in NeTCDF files. For the analyses that we conduct, there is little need to\nproduce outputs more frequently than that, bearing in mind that a single simulation\noutput at 50 year intervals is a NetCDF file of 80 GBs in size due to number of variables\nwe need to output and the large size of the grid. The ensemble of 100 simulations (plus\nprevious tests) thus represents a large amount of data to handle during post-processing\n(~10 TB).\n-Reviewer:\n“Line 169: (I would expand more on ensemble for clarification: is it ensemble of 100\nparameterized model simulations resulting from a group of parameter sets and /or a\ngroup of climate input used in the study)”\nAuthors:\nThank you for raising this point. The 10 chosen parameters lead to simulation-\ndependent modifications in both ice properties and boundary conditions (including\nclimate and SMB, but also basal conditions, bed topography, etc). More detailed\ninformation on this are provided in supplementary table 1, however we have added a few\ndetails to the relevant sentence for clarification:\n“Within this ensemble, we vary 10 key model parameters drawn from eight different\ncomponents of the glacier system which modify both ice properties and boundary\nconditions (see Supplementary Table 1). The 100 parameter sets are sampled between\ngiven parameter ranges using a Latin hypercube algorithm46, optimised with a maximin\ncriterion.”\n-Reviewer:\n“how is the initial range for each parameter selected? Is it from the output of the model\nsimulations or based on the model default values or from literatures?”\nAuthors:\nWe thank the reviewer for raising this important point.\nFor most parameters which are transferable between different models, the initial range\nof parameter values is selected to be a conservative range (in order to explore quite\nwidely the possible model responses) but bracketing values found to best represent the\nconditions of the AIF during the LGM, as found in the literature and by previous\nmodelling efforts (mostly from previous modelling work of Juergen Mey, Julien Seguinot\nand Guillaume Jouvet). For ensemble-varying parameters that are specific to IGM and\nthis experiment, however, the initial ranges were obtained by conducting a series of\nsensitivity tests before running the ensemble and bracketing reference values that\ngenerated a better model data fit.\nWe have added this information as a caption to Supplementary table 1 for increased\nclarity:\n“*The initial ranges of parameter values that are sampled from (by the Latin Hypercube\nalgorithm) are selected to be a conservative range (thus exploring quite widely the\npossible model responses) bracketing values found to best represent the conditions of\nthe AIF during the LGM, as found by previous modelling efforts, including the previous\nwork of Mey et al. (2016), Seguinot et al. (2018), Jouvet et al. (2023). For ensemble-\nvarying parameters that are specific to IGM and this experiment (i.e. the catchment-\nspecific precipitation offset parameter, and topographic control on yield stress\nparameters), initial ranges were obtained by conducting a series of sensitivity tests prior\nto running the ensemble and targeting ranges that bracket values shown to produce a\nbetter model data fit.”\n-Reviewer:\n“Line 225 -227: Any explanation for why the maximum AIF areas outside the range of\n156-165X 103 produce substantial under-or overestimations of the LGM ice extent\nreconstruction?”\nAuthors:\nAs part of our answer to the main comment of reviewer 1 on the model-data comparison\nscheme (now more quantitative) and the changes made to the relevant text (see\nseparate response), this statement is no longer part of the paper. However it is very\ninteresting to still address why this happens. We hypothesise this is likely related to LGM\nmoraines of the Alpine foreland being located on the edges of overdeepenings eroded by\nmultiple quaternary glaciations of similar cooling amplitudes, and which create sweet\nspots where modelled ice margins can easily stabilise due to reverse bed slopes and\nvoluminous moraine dams. However there likely exist threshold points whereby if the\nmodelled AIF is being pushed enough to grow more, the ice suddenly overspills these\nmoraine dams and is then able to flow much more extensively into the foreland due to\ndown-sloping beds. Contrastingly, if the model is not being forced to grow enough to\nthese sweet spots, then the reverse bed slopes cause most outlet glaciers to be far too\nretreated. The range of 156-165 x103 km2 is likely to represent the boundaries of these\nsweet spots: which are thus a consequence of the highly topographically-constrained\nnature of the AIF during the LGM.\n-Reviewer:\n“Line 228- 230: Could it be because of the deep learning algorithm that the model is\nproducing similar patterns for outlet glaciers despite extensive changes in the\nparameters? does this also mean these outlet glaciers control the evolution and\nbehaviour of the entire AIF?”\nAuthors:\nMany thanks for raising this really interesting point. As previously mentioned due to\nchanges made to the model-data comparison scheme (see separate reply to main\ncomments), this statement (which originally came from a semi-quantitative\nobservation) is not part of the paper anymore.\nHowever, we do not believe the iceflow emulator would be responsible for this, as it\nmodels ice flow with high fidelity to the Blatter-Pattyn solver (see Jouvet et al. 2023: IGM\nemulator paper) and is thus expected to produce very similar results to any traditional\nhigh-order model. Such a spatially-homogeneous response in modelled AIF-wide extent\nwas already observed in our previous modelling work (as part of the PISM modelling\nstudy of Jouvet et al. 2023), and we believe it is related to the highly topographically-\nconstrained nature of the AIF during the LGM. As the ice divides and ice flux patterns are\nnot able to change between simulations due to being heavily controlled by valleys and\ndrainage basins, the different extent and volume response of the AIF are heavily\ndominated by changes in the surface mass balance controlled by the PDD melt factors\nand refreezing parameters (see supplementary Table 1), whilst parameters controlling\nbasal conditions (thus sliding) and ice properties have relatively less impact. At the\nscale of the entire AIF (which is what is shown with these time series in Fig. 2), this leads\nto parallel time series and spatially-homogeneous responses in modelled AIF extent, as\nthose climate and PDD parameters are not space-dependent (the change from one\nsimulation to the next in “climate” or “SMB” forcing is the same across the domain). We\nalso agree that the evolution of main outlet glaciers and the rest of the AIF are thus\nhighly coupled.\n-Reviewer:\n“Line 241-243: how is this natural break for mean and standard deviation i.e. 180 m and\n100 m in elevation difference chosen? (is it from sieved simulation data set?)”\nAuthors:\nOriginally, it was found by observing sudden jumps in model-data fit in ice extent: where\ngroups of simulations more closely respected the empirical outline when under certain\nsieve values. However this was semi-quantitative as based on observations as well. We\nhave now changed the model-data comparison procedure to a more quantitative\napproach (see reply to main comment on this), and now choose the 3 sieve threshold\nvalues in order to incrementally remove 90% of the simulations and end up with a final\npool of 10 best simulations (10% of the ensemble) which we then conduct a final visual\ncheck on. This new approach results in the same 8 NROYs as were found before: as\nthese clearly stand out as the best-fit simulations and are indistinguishable in quality of\nmodel-data fit. This protocol has been updated in the text accordingly: in the results\nsection “model-data comparison scheme”:\n“Simulations with mean and standard deviation values exceeding 160 m and 100 m in\nelevation difference, respectively, were sieved out (n=32). This quantitative three-sieve\nmodel-data comparison procedure was designed to isolate a group of 10 best-fit\nsimulations (10% of the ensemble) for further inspection. Finally, a visual check of these\n10 simulations by three experts enabled to identify two outlier simulations which\nproduce far too extensive ice margins throughout the Western and Northwestern Alps.\nAfter removing these two outliers, the final eight Not-Ruled-Out-Yet simulations (NROYs)\nare found to be indistinguishable in quality and ice-extent fit with the empirical LGM\noutline and are thus used for all subsequent quantitative results reported in this study,\nusing NROY means and standard deviations in given metrics (Fig. 2).”\n-Reviewer:\n“Line 247: the best fit to moraines ensemble member selection: was how was this\nquantified? Visually?”\nAuthors:\nYes this is a great point. Simulation 37 can be isolated as our best-fit simulation almost\nentirely quantitatively. If increasing the harshness of our three sieves (i.e. in the new\nscheme: extent score > 0.84, mean of misfit with trimlines <148 m, stdev of misfit with\ntrimlines <96 m) incrementally and homogeneously across sieves such that the number\nof remaining simulations is progressively reduced: we basically end up with two last\nsimulations that present the best model-data fit scores: number 30 and 37. Number 30\nis one of the outliers which we identified during final visual check: due to far too\nextensive ice along the western and northern edges of the AIF.: and is thus clearly worse\nthan simulation 37 in model-data fit. Therefore, amongst the 8 NROYs identified (thus\nafter removing the two outliers), the amalgamated model-data fit statistics of simulation\n37 are better than the 7 others. However this simulation is only selected for visualization\npurposes: the 8 NROYs are pretty much indistinguishable in quality of model-data fit\n(when analysing them visually) and all quantitative results of the paper are thus reported\ntaking into account results from all 8 NROYs (more detail to the text has been added to\nclarify this, see response to reviewer 2’s comment on this topic).\n-Reviewer:\n“Line 280 – unprecedented is a very strong word, can it be justified? I would avoid using\nmore superfluous language personally.”\nAuthors:\nWe thank the reviewer for pointing this out and agree this sentence sounded a little too\ndramatic. We have now reworded to:\n“Overall, our model produces a thinner LGM ice field resulting in more valley-confined\nice flow, and leading to a substantial increase in model-data fit in ice surface elevation.”\n-Reviewer:\n“Line 326-327: why is that so? Is it because this outlet glacier is influenced by more local\nphenomena than the Alps-wide phenomenon and the ensemble parameter sets could\nnot capture the processes and drivers?”\nAuthors:\nWe thank the reviewer for this highly relevant question. Reviewer 4 asked a similar\nquestion related to what we think could explain the remaining model-data biases in both\nthickness and extent. We have provided three main groups of reasons for what we think\nto be the most impactful mechanisms in our reply and have added a paragraph to the\nresults section on this: see our reply to reviewer 4 further down. Also: here is the\nparagraph we added to the main text:\n“We believe remaining model-data misfits in ice thickness (+146 ± 12 m) and in the\nextent of certain outlet glaciers can either be, 1) Reduced further with higher-than-300 m\nresolution simulations and more extensive parameter-space explorations through larger\nensembles, 2) Related to unavoidable uncertainties inherent to paleo trimline and LGM\nmoraine identification due to challenges with landform preservation and dating, and 3)\nRelated to uncertainties and biases in our input climate and surface mass balance\nparameterisation. A combination of impact from these three mechanisms is also\nconsidered likely. For the latter, future improvements in the input climate through higher-\nresolution regional climate modelling and more transient simulations accompanied by a\nmore complex energy balance model may help reduce some of the remaining model-\ndata misfits.”\n-Reviewer:\n“Line 408: explain briefly what high-magnitude venturi effects mean (not necessarily all\nglaciologists know what it means)”\nAuthors:\nYes thank you for pointing this out, it is true that this term is mostly used in fluid\nmechanics and in climatology, not so often in glaciology (although very relevant, we\nbelieve). We have added this second sentence here:\n“When ice in such valleys drains large glacier catchments, NROYs produce greater flow\nconvergence, with higher-magnitude venturi effects causing ice flow speedups (Figs. 5,\n7). The venturi effect here applies as it describes the increase in an incompressible\nfluid’s velocity as it passes through a constriction, in order to respect the principle of\nmass conservation.”\n-Reviewer:\n“Line 557 – how good is the chronological evidence here?”\nAuthors:\nThe relevant sentence we wrote in the manuscript is: “The modelled Rhône outlet glacier\nreaches maximum extent around 24.6 ka, five centuries after most glaciers (including\nthe nearby Rhein outlet glacier) started to retreat from their maximum positions.”\nAnd here is a citation from the review by Ivy-Ochs (2015) about this specific margin:\n“10Be exposure dates for boulders of hornblende granite brought from the southern\nvalleys of the Valais located along the right lateral position indicate initial glacier\nwithdrawal no later than 24±2 ka (Ivy-Ochs et al., 2004)”\nNew geochronological constraints were also produced by Hofmann et al. (2024), Graf et\nal. (2015), and Rey et al. (2017 and 2020).\nHere is a citation from Hofmann et al. (2024) regarding outlet glaciers of this region\n(Rhone, Reuss):\n“According to optically stimulated and infrared-stimulated luminescence ages, glacio-\nfluvial and glacio-lacustrine sediments near the MIS 2 maximum position of the Reuss\npiedmont glacier (Fig. 5) were deposited at 25.1 ± 2.4 ka and 24.2 ± 2.2 ka, respectively\n(Gaar et al., 2019). Glacier retreat from the morphologically indistinct moraines at the\nmaximum MIS 2 position of the Reuss glacier, the Untertannwald ice-marginal position,\nwas underway by 21.3 ± 1.0 ka”\nThus, the fit in timing with our modelled LGM in this region seems ok. However please\nnote that our model-data comparison scheme is time-independent (as shown on Fig.\n1a) and only assesses the fit in LGM AIF geometry with moraines generally described by\nthe community as “LGM in timing”, but that could have been deposited anytime\nbetween our simulation’s start and end times (35 to 18 kyr BP): which safely brackets the\nentire LGM period. As this is not our objective, in this work, we do not assess the fit in the\nprecise timing of our modelled LGM with dates attributed to LGM margins in each region.\nWe have included more details relevant to this in our reply to the main comment by\nreviewer 1 on model-data comparison: and have also made sure to add this to the\nmodel-data comparison scheme for more clarity:\n“Here, the modelled LGM state for each simulation is obtained by computing the\nmaximum thickness reached in each pixel at any time during the simulation, thus\nproviding a time-independent map of maximum modelled thickness and extent (Fig. 1a).\nThe precise dated timing of the LGM in specific regions is thus not directly compared\nwith the timing of the modelled LGM as part of this model-data comparison scheme,\nwhich exclusively aims at finding simulations with the most data-consistent AIF\ngeometry (see Supplementary Fig. 19).”\n-Reviewer:\n“Line 637 – how else could it inform future investigations – see comments above.”\nAuthors:\nWe thank the reviewer for this comment. Please see our reply to the main reviewer 1\ncomment on the wider implications section of the discussion: to which we have now\nmodified and added more text in line with the reviewer’s suggestions.\n-Reviewer:\n“Line 647 – there are also implications for data as well”\nAuthors:\nWe strongly agree with this, please see our reply to the main reviewer 1 comment on the\nwider implications section of the discussion.\n-Reviewer:\n“Line 656 – this is very vague and I suggest rewriting this paragraph. Greenland &\nWAIS/EAIS have very few nunataks and requires a different approach. Can a very small\nnumber of spatially limited dispersed nunataks improve modelling here?”\nAuthors:\nWe apologise that these few sentences were not very clear. Originally, we were here\nreferring much more broadly (moving away from the Alps and model-data comparison\nwith trimlines, nunataks etc) to the new technology of GPU-based glacier evolution\nmodels and physics-informed machine learning (which is here applied for the first time\nto a continental scale ice body), and which permits significant gains in computational\nmeans and thus in model resolution. What we referred to in this sentence is that\nmodelling the past and future of the Greenland Ice Sheets and Antarctic Ice Sheets,\nbesides having huge future implications for society, is currently limited by\ncomputational costs and thus model resolutions (especially in paleo timescale). The\nAntarctic Ice Sheet, for instance, in paleo timescales is usually modelled at 15 km\nresolution at best (in the interior). At the 15 km, the interactions between the ice sheet\nand its bed are very likely to be mis-represented. We know that for these ice sheets,\ncalibrating the models over paleo millennial timescales is important and impacts future\nprojections. Further developing GPU-based models may bring significant improvements\nin that regard, enabling to model the entire ice sheets at sub-km resolution with\nhomogeneous grids (we anticipate this will soon be possible thanks to the new\ngeneration of GPUs).\nTo increase the clarity of that message, we have reworded and increased the length of\nthat particular paragraph to:\n“We also believe this new technology may permit more accurate simulations of past and\nfuture Greenland- and Antarctic-ice-sheet change, by facilitating ice-sheet-wide\nsimulations at considerably higher resolutions than previously achieved and at lower\ncomputational costs. This would allow for broader explorations of parameter spaces and\nice-sheet evolution scenarios, with implications for more accurately modelling basal\nconditions and ice-ocean interactions towards grounding lines, crucial for better\nprojecting the future response of continental ice sheets to climate change38. Moreover,\nreducing Quaternary model-data misfits, as is achieved here with the AIF, may be\nimportant for improving the paleo initialisation procedures of Antarctic and Greenland\nice sheet models, which greatly impact their future projections38,41. Therefore, we\nexpect physics-informed, GPU-based models to bring advances in projecting future\nglobal sea level rise and environmental change associated with the response of\ncontemporary ice sheets to climate change.”\nPlease also refer to our reply to the main reviewer 1 comment on the wider implications’\npart of the discussion.\n-Reviewer:\n“Line 659 – a rather superficial sentence to end on. You can finish with a more exciting\ncall to arms and set the agenda for the community here.”\nAuthors:\nWe thank the reviewer for this comment and have added a new conclusion sentence\nthat attempts to do more agenda-setting for the community:\n“To conclude, we believe the development and application of physics-driven AI and\nGPU-based glacier and ice-sheet models should become a primary avenue of future\nresearch due to its wide and multidisciplinary implications for the fields of Quaternary\nand cryosphere science.”\n######################################################################\nReviewer 2\n-Reviewer:\n“Abstract (no line numbers) just some word-smithing suggestions, results ‘have’\nimplications instead of ‘yield’ implications. Suggest saying something about how the\ndecreased computation time allows users to model glaciers with ensembles, the way to\ngo to test parametrizations, allows for uncertainty quantification, the list goes on. This\ngoes well beyond modeling glaciers in the Alps better.”\nAuthors:\nWe thank the reviewer for this point. We have rephrased to ‘have” implications, and have\nadded some more details to the last sentence of the abstract:\n“Furthermore, this study demonstrates that physics-informed AI-driven glacier evolution\nmodels can overcome the bottleneck of high-resolution continental-scale modelling,\nand enable to better explore parameter spaces through reduced computational costs,\nboth of which are required by glacier and ice-sheet modelling to accurately describe\ncomplex topographies, ice dynamics, and test model parameterisations.”\n-Reviewer:\n“Line 7, ‘one of the better reconstructed’ reword..”\nAuthors:\nWe rephrased this to: “making the extent of the AIF during the Last Glacial Maximum\n(LGM: ~30-19 ka5,6) one of the most well known in the world”\n-Reviewer:\n“Line 12, ‘also help reconstruct’ reword…”\nAuthors:\nThis was rephrased to: “has implications for understanding such mechanisms and can\nalso improve reconstructions of late-Quaternary human and ecological history”.\n-Reviewer:\n“Line 28 – these trimlines are debated all over the world (not just Britain and\nScandinavia, but of course you can’t cite all)”\nAuthors:\nWe thank the reviewer for this point, which we fully agree with. Of course we are limited\nin the number of reference which we have already reached the limit of. But here we have\nrephrased to clearly state that Scandinavia and Britain are only two examples:\n“as was previously debated in other glaciated regions, such as Britain and Scandinavia\nfor instance23.”\n-Reviewer:\n“Line 31 – this paragraph sure has a lot of acronyms. Acronyms are overrated (!), do we\nreally need them at all anymore? Are there space constraints and word counts? Suggest\ndoing away with them, it just makes writing more inaccessible and difficult to\nunderstand.”\nAuthors:\nWe thank the reviewer and we fully agree with this comment, it would be nice to move\naway from acronyms as much as possible. However here we are also limited by word\ncount: with a maximum of 5000 words for the main text and 3000 words for the Methods.\nAs of now, the main text is at 4800 words, so a bit tight already. The IGM, AIF, and LGM\nacronyms are used more than 50 times each across the manuscript and would thus\nincrease the word count far beyond the limit if spelled out.\nThe acronym GPU is actually very important to keep (unfortunately). It is the key to the\nentire new modelling approach followed in this study. Being able to compute on the GPU\nis what makes the computational gains possible. We believe this term will become\nincreasingly important within environmental sciences in general and people need to\ntherefore know what it means. Unfortunately, computer scientists only refer to GPUs\nusing the acronym, and never the full term “Graphics processing unit”.\nWe have however removed the acronym ‘CPU’ from the paper, which was not used\nenough to be relevant.\n-Reviewer:\nLine 47, this reference of 1-2 km grid should come several sentences earlier when\ndiscussing that older model\nAuthors:\nThere is in fact a reference to this resolution earlier in the introduction (originally line 23):\nwith this sentence:\n“This ice thickness overestimation has been attributed to oversimplified ice flow\ndynamics and the coarse spatial resolutions of previous AIF-wide models,\ncomputationally restricted to 1-2 km”.\nHowever, the term “previous” was not there and we here added it following the\nreviewer’s comment to make it clearer that these were indeed the best-achieved\nresolutions for previous models of the Alpine Ice Field before our study.\n-Reviewer:\nLine 58, upon reflection, half or more than half of the introduction is just a summary of\nthe results and the paper – that’s not really introduction. Yes, many times introductions\nclose with perhaps a final paragraph on ‘what we did’, but this seems a bit much.\nSuggest reducing summary text that appears in the introduction, and in place adding\nmore actual intro material. As is, it’s like reading results, then going to the results section\nto read again the same results..\nAuthors:\nWe thank the reviewer for this comment. Indeed the introduction is structured in three\nparagraphs. The first two describe the context of AIF reconstruction, and the issues with\nprevious modelling efforts and the model-data misfits, respectively (thus highlighting\nthe context and motivation for our work). The third and last introduction paragraph,\nhowever, was written as a summary of the new approach followed in our work, and the\nkey results it produces. We agree with the reviewer this third paragraph could be made\nshorter and include less results, as it needs to make the reader want to read the rest of\nthe paper without giving away too much detailed information.\nWe have therefore removed a few sentences from the introduction paragraph 3 regarding\nspecific results. In the new text, results are now only summarized in two sentences (with\nno specific numbers given, only trends). The rest of the third paragraph is more focused\non introducing the nature of the novel method employed (GPU-based modelling) and its\nwider implications. The latter has been expanded a bit in line with new changes to the\ndiscussion following comments by reviewer 1 :\n“To address this conundrum, we apply a novel approach by modelling the AIF with the\nInstructed Glacier Model (IGM)31,32, a thermo-mechanically coupled three-dimensional\nglacier evolution model. IGM makes use of recent improvements in physics-informed\nmachine learning to accelerate the high-order “Blatter-Pattyn” solver of ice flow33,\nenabling efficient computation on graphics processing units (GPU). The resulting model\nsimulates glacier motion with fidelity while offering the same ability to perturbate\nparameters controlling ice properties as traditional central-processing-unit-based\nmodels, but at a fraction of the computational cost. Here, we use IGM to simulate the\nAIF’s transient evolution from 35 to 18 ka, thus bracketing the full LGM period, at a\nspatial resolution of 300 m, an order of magnitude higher than previously achieved20.\nComputationally unfeasible with traditional glacier models, our perturbed parameter\nensemble experiment shows that increasing spatial resolution substantially reduces\nmodel-data misfits in both LGM ice extent and thickness across the Alps, resulting in a\nthinner ice field than previous 1-2 km resolution models20,22. This new approach enables\nus, for the first time, to produce a high-resolution, data-consistent, yet physics-based\ntransient simulation of the AIF’s last glaciation. The results have wider implications for\naddressing open research questions on Quaternary glacial erosion processes,\npostglacial isostatic rebound, and paleoclimate in the European Alps. Our study also\nreveals that GPU-based glacier evolution modelling is a promising tool for better\nreconstructing other paleo ice fields, but also for past and future ice-sheet-scale\nmodelling. Indeed, by substantially reducing computational costs, this new approach\npermits advances in modelling resolution and parameter space exploration, both\nessential to resolve complex topographies and model more accurate past and future\nglaciated environments.”\n-Reviewer:\n“Line 107, not sure I understand the “Are shown; blah” syntax… Change?”\nAuthors:\nMany thanks for pointing this out, we have now changed the wording to:\n“Panel a shows the 2-D field of modelled time-independent maximum ice thickness for\nsimulation 37”\n-Reviewer:\n“Line 143, are you sure “results” is the best description of these first paragraphs – seems\nmore like methods to me. Can you come up with a title that’s neither, like “the Instructed\nGlacier Model” or something?”\nAuthors:\nWe here provide the same reply as for a comment by reviewer 4 on a very similar point.\nYes we agree this can be a bit difficult. As far as we understand the main text of the\npaper should give all the tools for the readers to understand what the experiment is, how\nit was conducted, the results it produced and their wider implications. More precise\ndetails on the detailed techniques employed should be moved to the methods.\nHowever, some methods are simply too unavoidable to not be mentioned in the main\ntext, otherwise the experiment cannot be understood. This is the case with this study\nwhich relies on a new approach and new model. This means that parts of the ‘results’\nend up being more methods-like. This is why we used explicit headings in this case: with\nparagraphs being called “High resolution model setup” and “Model-data comparison\nscheme”. This slight difference with more conventional papers occurs in most Nature\nCommunications papers, we find. The editorial guidelines on how to name different\n“main” paper sections is fairly strict, it seems, and we don’t believe it is possible to\nmove text between the “introduction” and “results” sections in Nature Communication\npapers, unless the editor says otherwise and allows it to occur in this case?\nWe thus let the editor make a decision on this, based on whether it is possible to change\nthis in Nature Communications papers or not.\n-Reviewer:\n“Line 149, ‘follows and improves’ well I hope it would improve! It sort of goes without\nsaying, can you replace with ‘builds on’ ?”\nAuthors:\nWe agree that “builds on” would be more appropriate and have now reworded the\nsentence using that instead, many thanks for pointing this out.:\n“we design a high-resolution (300 m) model setup with IGM24,25 that builds on the Parallel\nIce Sheet Model (PISM27) experimental setup of Jouvet et al.10.”\n-Reviewer:\n“Line 156, need word ‘materials’?”\nAuthors:\nIndeed this word is not necessary and has now been removed from the sentence.\n-Reviewer:\n“Line 166, In ‘the’ main overdeepened valleys…”\nAuthors:\nThanks for spotting this, we added the word accordingly.\n-Reviewer:\n“Line 190 – figure 2/study design… I recommend explaining in the text where the sieve\nboundaries come from. How are they derived? Why not choose a narrower range? Or a\nbroader range? Is it arbitrary, or based on something? Sorry if I missed it.”\nAuthors:\nYes this is a very good and important point. Thank you. Please refer to our longer reply to\nreviewer 1 on the model-data comparison scheme, which has now been modified to\nbecome more quantitative and robust, we believe. This new scheme moreover produces\nthe same set of 8 NROYs which are found to be indistinguishable in model-data fit\nquality, which is re-assuring as it shows that slightly different methods in model-data\ntests (and different sieves) lead to similar results. Increasing the harshness of the sieves\nfurther would result in the removal of simulations which are (visually) indistinguishable\nin their fit with the empirical LGM outline, while their fit with the trimlines are also very\nsimilar: and would thus be difficult to justify.\n-Reviewer:\n“Line 203, can you add a comma in that big number (6 digits)”\nAuthors:\nWe have changed this accordingly.\n-Reviewer:\n“Line 207, can you do away with the footnote and just add the sentence to the\nparagraph?”\nAuthors:\nWe thank the reviewer for pointing this out and have now added the sentence and\nremoved the footnote accordingly:\n“While running a 20 kyr-long simulation at 300 m would require ~2.5 years with PISM\nusing a 32-core (3.7 GHz) computer (or ~6 months using a state-of-the-art high-\nperformance computer), it is achieved in only ~2.5 days with IGM using a single GPU\n(Nvidia RTX 4090).”\n-Reviewer:\n“Line 212, recommend instead of ‘we run a series’ write ‘we ran 100 simulations from 35\nto 18 ka’. I think at this point we know it’s using the 300 m grid.”\nAuthors:\nThis sentence was changed accordingly.\n-Reviewer:\n“Line 213, data-model agreement is mentioned, but I, for one, didn’t understand yet\nwhat data – were they introduced yet, or mentioned yet?”\nAuthors:\nMany thanks for pointing this out: the model-data comparison aspect is presented for\nthe first time in the next paragraph, which follows 3 sentences later. We thus added “see\nbelow” behind “model-data agreement” in a bracket within this sentence.\n-Reviewer:\n“Line 226 – how is this range chosen?”\nAuthors:\nThank for for raising this very good point. However due to the new change made to the\nmodel-data comparison scheme (see detailed reply to reviewer 1’s main comment),\nwith a first sieve modified to be a fully quantitative assessment and scoring of model-\ndata fit in LGM ice extent, we now do not have to specify a range for this anymore, and\nthis sentence has now been removed.\n-Reviewer:\n“Line 233, might be nice here to add that 18 sims are left”\nAuthors:\nIn the new text that describes the new model-data comparison scheme, we have made\nsure to make this more explicit indeed:\n“This quantitative three-sieve model-data comparison procedure was designed to\nisolate a group of 10 best-fit simulations (10% of the ensemble) for further inspection.”\n-Reviewer:\n“Line 244, say ‘3-step’ but that confused me because text describes it as a 2-step. I\nrealize there is a volume step that is derived from ice-thickness (?) so it’s kind of like\nthere are two parts of a second step – but I found the description of this as a ‘3-step’\nprocess a little confusing. Maybe wording will help.”\nAuthors:\nWe agree this was worded in a confusing way before. In the new relevant sentence\nshown directly above in reply to the previous comment: please note that we have\nreworded using “three-sieve” rather than “three-step”, which is more explicit.\n-Reviewer:\n“Line 246, simulation 37. I agree with your choice to pick one for illustration. I got hung\nup on this a little bit, because if, in the end, you just pick one among the 100 that\nmatches the data best, then why do the 3-step sieve process? Why have 8 NROYs and\nnot just pick the best? I think you have 8 NROYs (not just one) so you can have some\nattempt at assigning an uncertainty to your result, no? Maybe you can articulate this\nbetter, if that’s the case. Otherwise, ignore the sieving and just choose the single best.”\nAuthors:\nYes indeed we use the 8 indistinguishable NROYs to use for better uncertainty\nquantification and to find potential patterns (or not) in most suitable parameter values,\nbut also to not be biased towards one parameter configuration. This is why, in the\n“Improved LGM model-data fit” and throughout the discussion section, we\nsystematically report our results using NROY-mean and NROY-stdev statistics (with ctrl\n+ f you will find the term “NROY mean” comes up more than 10 times throughout the text\nand captions). However we agree this was not made clear enough in this “Model-data\ncomparison scheme” section.\nFor more clarity, we have thus added the sentence:\n“the final eight Not-Ruled-Out-Yet simulations (NROYs) are found to be indistinguishable\nin quality and ice-extent fit with the empirical LGM outline and are thus used for all\nsubsequent quantitative results reported in this study, using NROY means and standard\ndeviations in given metrics”\n-Reviewer:\n“Line 274, ‘thus permitting to compute’ reword..”\nAuthors:\nWe have rephrased to: “Note panel a only displays coloured pixels where ice is modelled\nby both IGM and PISM simulations (when maximum ice thickness is reached in each\npixel), as this is required to compute the thickness difference”.\n-Reviewer:\n“Line 283, reword for clarity”\nAuthors:\nWe agree this was a little confusing. We have reworded this sentence to:\n“The positive bias in modelled ice thickness is reduced by 200% and 450% compared\nwith these earlier studies, respectively.”\n-Reviewer:\n“Line 293, suggest phrasing ‘… increasing model resolution from x km to x m\nsystematically …”\nAuthors:\nThis has now been changed accordingly.\n-Reviewer:\nLine 297, resulting ‘thickness difference’ ?\nAuthors:\nWe have now replaced “thinning pattern” by “thickness difference” as suggested.\n-Reviewer:\nLine 302, replace ‘modelled’ with ‘achieved’\nAuthors:\nWe have now replaced it accordingly, many thanks for pointing this out.\n-Reviewer:\nLine 308, remove comma\nAuthors:\nNicely spotted, many thanks. Now removed.\n-Reviewer:\nLine 308/309, huh? Reword?\nAuthors:\nYes indeed this sentence was a little confusing. We have reworded to:\n“Across our 100 simulations, we find highly variable modelled basal ice temperatures at\nthe locations of trimlines during the LGM, with no apparent clusters towards distinct\ntemperatures. This suggests no relationship exists in our modelling between thermal\nboundaries of basal ice and trimline formation (see Supplementary Fig. 4).”\n-Reviewer:\n“Line 329, However, we…”\nAuthors:\nChange made accordingly.\n-Reviewer:\n“Line 331, ‘due to too extensive…’ reword”\nAuthors:\nWe have rephrased to: “However, we note three instances, i.e. the Rhône, Dora-Baltea\nand Isère outlet glaciers, where our NROYs model too extensive ice margins leading to a\nworse fit with data than Jouvet et al10’s model (Fig. 4a)”\n-Reviewer:\n“Line 402, are these two phenomena not just one?”\nAuthors:\nMany thanks, we now realise this sentence was not 100% clear. By “these two\nphenomena”, we meant the higher resolution modelling and the reduction in ice\nthickness. For more clarity, we rephrased to:\n“Our higher-resolution model simulations (300 m instead of 1-2 km in previous\nstudies8,10) result in thinner ice and glacier surface lowering across the AIF during the\nLGM. We hypothesise that localised ice flow speedup is the main underlying mechanism\nresponsible for this result.”\n-Reviewer:\n“Line 470, suggest ‘cold based regions cover 14±3% of the AIF (NROY mean)(Fig. 6b),\nmore than the 11% modelled at 2 km resolution’”\nAuthors:\nWe have made this change accordingly.\n-Reviewer:\n“Line 581, change ‘positive gradients’ to ‘increases’”\nAuthors:\nWe have made this change accordingly.\n-Reviewer:\n“Lines 695 and 890, find a better word/phrase to replace ‘hardly’ with, which is a bit\ncolloquial.”\nAuthors:\nMany thanks for spotting this,\nWe have now reworded to : “with negligible accuracy losses”\nAnd to: “The impact of removing these 43 presumed outliers on the mean model-data\nmisfit between modelled ice surface elevations and reported trimline elevations is\nnegligible (on average a ~2% reduction)”\n######################################################################\nReviewer 4\n-Reviewer:\n“L34: \"so-called high-order (three-dimensional)\". I'm not sure that any of these terms are\nrelevant to a general audience - these terms (or the subtexts associated with them) are\nterms-of-art for technical glacier modelers and may not need to appear here.”\nAuthors:\nWe agree that the term “high-order” is most likely too technical for a general audience\nand thus removed it from the paper Abstract. In the introduction, we also simplified the\nsentence “to accelerate the so-called high-order (three-dimensional) “Blatter-Pattyn”\nsolver of ice flow”, to replace it by “to accelerate the high-order “Blatter-Pattyn” solver of\nice flow”. We however consider it important to mention that IGM is informed by high-\norder physics, and not the SIA or a hybrid (SIA/SSA) scheme, as this matters for\nobtaining better results in complex topographies, as mentioned in the results and\ndiscussion sections. Moreover, we believe most glaciologists (including non-modelers)\nare now exposed enough to modelling studies to know the difference between SIA, high-\norder, and full-stokes ice flow physics.\n-Reviewer:\n“L37: What does \"identical levels of parameter control\" mean?”\nAuthors:\nWe thank the reviewer for pointing this out, this sentence is indeed not very clear. Here\nwe simply mean to say that IGM, despite being a slightly different model than traditional\nnumerical glacier models (CPU-based) due to the ice flow solver being emulated by a\nneural network, still provides the same ability to play with the same parameters of ice\ndynamics and ice flow motion as other models. In other terms, one can tweak exactly\nthe same knobs as before.\nTo make this clearer, we modified the relevant sentence to:\n“The resulting model simulates glacier motion with fidelity while offering the same ability\nto modify parameters controlling ice properties, but at a fraction of the computational\ncost of traditional central processing unit (CPU)-based models.”\n-Reviewer:\n“L159: Awkward sentence.”\nAuthors:\nWe thank the reviewer for pointing this out. For more clarity, we reworded and modified\nthe sentence to:\n“The experiment setup of Jouvet et al.10 is however improved by using an Alps-specific\nclimate record (instead of the EPICA34 ice core record) as input to our glacial index\nscheme10,35 (Fig. 2e). The applied signal combines the Bergsee lacustrine record36 (35-30\nka: Black Forest, Southern Germany) and the Sieben Hengste speleothem δ18O record37\n(30-18 ka: Bernese Alps, Switzerland).”\n-Reviewer:\n“L146--247: I appreciate that Nature-type journals don't want methods in the\nmanuscript, which is why these sections (which clearly are methods) are labeled as\n`results'. But maybe it's okay to just call them methods anyways? Or have an\n`experimental setup' section before results?”\nAuthors:\nYes we agree this can be a bit difficult. As far as we understand the main text of the\npaper should give all the tools for the readers to understand what the experiment is, how\nit was conducted, the results it produced and their wider implications. More precise\ndetails on the detailed techniques employed should be moved to the methods.\nHowever, some methods are simply too unavoidable to not be mentioned in the main\ntext, otherwise the experiment cannot be understood. This is the case with this study\nwhich relies on a new approach and new model. This means that parts of the ‘results’\nend up being more methods-like. This is why we used explicit headings in this case: with\nparagraphs being called “High resolution model setup” and “Model-data comparison\nscheme”. This slight difference with more conventional papers occurs in most Nature\nCommunications papers, we find. The editorial guidelines on how to name different\n“main” paper sections is fairly strict, it seems, and we don’t believe it is possible to\nmove text between the “Introduction” and “Results” sections in Nature Communication\npapers, unless the editor says otherwise and allows it to occur in this case?\n-Reviewer:\n“L241: Slightly more detail on what constitutes a `natural break' and how robust that\nselection is would be helpful for reproducibility.”\nAuthors:\nWe thank the reviewer for pointing this out: this section of text has now been updated\nand modified as the model-data comparison and sieving scheme has changed, in\nresponse to reviewer 1’s main comment (see more detailed reply there). We do not aim\nto identify these natural breaks anymore but instead aim to isolate NROYs that are\nindistinguishable in quality of model-data ice extent fit with the empirical LGM outline.\nThese natural breaks are not mentioned anymore in the text. The new scheme is more\nquantitative and thus more reproducible than the original approach, we believe. (see our\nreply to main comment by reviewer 1 on this).\n-Reviewer:\n“L279: What does `topographically constrained' mean?”\nAuthors:\nWith this term we here meant that the modelled AIF, because thinner than before, is now\nmore coupled to the topography (glacier catchments = fluvial catchments) thus\nbehaving like an ice field as opposed to an ice sheet that is thick enough to decouple\nfrom the topography. Indeed, this may not be well understood if termed “topographically\nconstrained”. We thus reworded this sentence to:\n“Overall, our model produces a thinner LGM ice field resulting in more valley-confined\nice flow, and leading to a substantial increase in model-data fit in ice surface elevation.”\n-Reviewer:\n“L282: It seems that the model remains considerably biased relative to observations,\nimprovement versus previous works notwithstanding. I think it would be of value to try to\nexplain why this remaining bias exists - especially relative to mean ice thickness or\nsomething - and if it's something that needs to be considered when contextualizing\nresults.”\nAuthors:\nIndeed some (although small relative to previous studies) remaining bias in the extent of\ncertain outlet glaciers and of +146 m in ice thickness between our best-fit models and\nempirical reconstructions can be observed. We believe it can mostly be attributed to\nthree possible mechanisms, or a combination of them.\n-1) the remaining ice thickness bias will likely be further reduced if running the\nsimulations at even higher spatial resolution than 300 m and if running larger ensembles\nenabling a deeper parameter space exploration. At the moment we are unable to run at\nhigher resolution due to limited GPU memory, but this will shortly be made possible by\nnew generations of GPUs (such as the NVIDIA H100) which will undoubtedly permit to\nmodel the entire AIF at 100 m and possibly to run larger ensembles. We don’t expect this\nthickness offset to disappear, but we don’t know until we can try. The full complexity of\nthe bed topography, especially over steep slopes, is still not perfectly well represented\nusing a 300 m resolution grid. Significantly more realistic ice flow, but also avalanche\nand snow redistribution processes over steep slopes and rock cliffs would be obtained\nwith a 100 or 50 m topography, which may help reduce the ice thickness offset further in\nsome cases.\n-2) Observing trimlines from bedrock surfaces transitioning from ice moulded erosion\nregimes to weathered and frost-shattered erosion regimes can be challenging and yields\nuncertainties, especially due to post deglaciation processes that can disturb, erode, and\nremove the highest elevated ice-moulded bedrock surfaces. The fact that the remaining\nmisfit is only of +~150 m and the R2 of ~0.92 (Fig. 1b) strongly suggests the assumption\nthat trimlines are formed near the ice surface is valid, in the European Alps at least.\nHowever it does not mean a trimline observation should necessarily be considered the\nexact ice surface. What we mean here is, not all uncertainties are necessarily linked to\nthe modelling , and it is possible a remaining misfit of ~100 m is unavoidable due to\nnature of the trimline observation data (for paleo reconstructions) and may represent an\nacceptable range of uncertainties in this type of model-data comparison. This is just a\nhypothesis, however, time and further investigations will help validate or refute this.\n3) Although the input climate used in this study is likely the best available input climate\ndata for the LGM in the European Alps (Russo et al., 2024), it still only comes from a\nsingle regional climate model (WRF), and is just one snapshot in time (the LGM at 24 kyr\nBP or the PI), and is still 2 km resolution, which may not be enough to capture the full\ninter-valley variations in past climate and mountain-specific climatology. Future\nimprovements in the input climate through more regional climate modelling\nreconstructions with more transient climate simulations and at higher spatial resolution\nwill most likely help reduce some of the remaining biases in ice thickness. Once the\nlargest uncertainties linked to regional climate modelling of the LGM are reduced, a\nbetter SMB model than the PDD (e.g. energy balance models) could also help better\nrepresent some of the missing surface ablation mechanisms (e.g. sublimation).\nBased on these thoughts, we have added a new paragraph to the relevant section of the\nresults (end of “Improved LGM model-data fit” section) to contextualize the results a bit\nmore:\n“We believe remaining model-data misfits in ice thickness (+146 ± 12 m) and in the\nextent of certain outlet glaciers can either be, 1) Reduced further with higher-than-300 m\nresolution simulations and more extensive parameter-space explorations through larger\nensembles, 2) Related to unavoidable uncertainties inherent to paleo trimline and LGM\nmoraine identification due to challenges with landform preservation and dating, and 3)\nRelated to uncertainties and biases in our input climate and surface mass balance\nparameterisation. A combination of impact from these three mechanisms is also\nconsidered likely. For the latter, future improvements in the input climate through higher-\nresolution regional climate modelling and more transient simulations accompanied by a\nmore complex energy balance model may help reduce some of the remaining model-\ndata misfits.”\n-Reviewer:\n“L408: Why would it be the case that the inclusion of longitudinal stresses produces\nhigher velocities?”\nAuthors:\nYes indeed we see where the confusion could lie here. In this sentence, we meant more\nthat due to Glen’s flow law, a high-order scheme such as the “Blatter-Pattyn” model\nwhich includes three-dimensional coupling of both the vertical and horizontal\ncomponents of strain will necessarily lead to a lower ice viscosity under a given stress\nthan when only considering the vertical (SIA) or horizontal (SSA) components of strain\nseparately.\nWe have rephrased this to:\n“Here, the high-order Blatter-Pattyn solver26 used in IGM25 considers both the vertical\nand horizontal components of strain through the ice column49, which helps produce\nstrong ice velocity gradients and promotes the formation of lateral shear margins in\nvalley-confined settings, thus encouraging higher glacier velocities when increasing\nspatial resolution50.”\n-Reviewer:\n“L676: What is the numerical method used to update the ice geometry? My\nunderstanding is that IGM only emulates the velocity field. If that's not the case, then\nhow is it ensured that essential conservation laws (especially mass) are upheld?”\nAuthors:\nWe thank the reviewer for pointing this out. It seems we forgot to write a key sentence in\nthis part of the method which explains this. It is indeed true to say that IGM only\nemulates the ice velocity field.\nWe’ve now added this important sentence to this paragraph:\n“Equation (1) is solved using an explicit upwind finite volume scheme on a regular grid\nwhich allows the model to update the ice thickness while conserving mass.”\n######################################################################"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Review of Leger et al., R1\nMany thanks for the thorough reubuttal and the swift resubmission of the manuscript. I have read through the authors\nresponse and the revised manuscript and can see the hard work that has gone into doing this.\nIn response to our queries the authors have proposed a new method of comparing the model and the data, which is\nappropriate and an improvement on the original method.\nThe clarifications requested have all been attended to. There is the requested added detail in the supplementary information\nand the wider discussion improved. Each comment has had a careful and thoughtful response.\nI have no further comments for the improvement of this paper and look forward to seeing it published.", "reviewer_2": "The authors did a fine job responding to my comments. I have nothing further.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_4": "Thank you for the carefully considered responses to my previous comments. I have no further suggestions."}, "author_response": ""}]

Decision letter