Our limited knowledge of the climate prevailing over Europe during former glaciations is the main obstacle to reconstruct the past evolution of the ice coverage over the Alps by numerical modelling. To address this challenge, we perform a two-step modelling approach: First, a regional climate model is used to downscale the time slice simulations of a global earth system model in high resolution, leading to climate snapshots during the Last Glacial Maximum (LGM) and the Marine Isotope Stage 4 (MIS4). Second, we combine these snapshots and a climate signal proxy to build a transient climate over the last glacial period and force the Parallel Ice Sheet Model to simulate the dynamical evolution of glaciers in the Alps. The results show that the extent of modelled glaciers during the LGM agrees with several independent key geological imprints, including moraine-based maximal reconstructed glacial extents, known ice transfluences and trajectories of erratic boulders of known origin and deposition. Our results highlight the benefit of multiphysical coupled climate and glacier transient modelling over simpler approaches to help reconstruct paleo glacier fluctuations in agreement with traces they have left on the landscape.

Convolutional neural networks (CNN) trained from high-order ice-flow model realisations have proven to be outstanding emulators in terms of fidelity and computational performance. However, the dependence on an ensemble of realisations of an instructor model renders this strategy difficult to generalise to a variety of ice-flow regimes found in the nature. To overcome this issue, we adopt the approach of physics-informed deep learning, which fuses traditional numerical solutions by finite differences/elements and deep-learning approaches. Here, we train a CNN to minimise the energy associated with high-order ice-flow equations within the time iterations of a glacier evolution model. As a result, our emulator is a promising alternative to traditional solvers thanks to its high computational efficiency (especially on GPU), its high fidelity to the original model, its simplified training (without requiring any data), its capability to handle a variety of ice-flow regimes and memorise previous solutions, and its relatively simple implementation. Embedded into the ‘Instructed Glacier Model’ (IGM) framework, the potential of the emulator is illustrated with three applications including a large-scale high-resolution (2400x4000) forward glacier evolution model, an inverse modelling case for data assimilation, and an ice shelf.

Future mass loss projections of the Greenland ice sheet require understanding of the processes at a glacier terminus, especially of iceberg calving. We present detailed and high-rate terrestrial radar interferometer observations of Eqip Sermia and Bowdoin Glacier, two outlet glaciers in Greenland with comparable dimensions and investigate iceberg calving, surface elevation, velocity, strain rates and their links to air temperature, tides and topography. The results reveal that the two glaciers exhibit very different flow and calving behaviour on different timescales. Ice flow driven by a steep surface slope with several topographic steps leads to high velocities, areas of extension and intense crevassing, which triggers frequent but small calving events independent of local velocity gradients. In contrast, ice flow under smooth surface slopes leaves the ice relatively intact, such that sporadic large-scale calving events dominate, which initiate in areas with high shearing. Flow acceleration caused by enhanced meltwater input and tidal velocity variations were observed for terminus sections close to floatation. Firmly grounded terminus sections showed no tidal signal and a weak short-term reaction to air temperature. These results demonstrate reaction timescales to external forcings from hours to months, which are, however, strongly dependent on local terminus geometry.

Data assimilation in high-order ice flow modeling is a challenging and computationally costly task, yet crucial to find ice thickness and ice flow parameter distributions that are consistent with ice flow mechanics and mass balance while best matching observations. Failing to find these distributions that are required as initial conditions leads to a disequilibrium between mass balance and ice flow, resulting in nonphysical transient effects in the prognostic model. Here we tackle this problem by inverting an emulator of the Stokes ice flow model based on deep learning. By substituting the ice flow equations using a convolutional neural network emulator, we simplify, make more robust and dramatically speed up the solving of the underlying optimization problem thanks to automatic differentiation, stochastic gradient methods and implementation of graphics processing unit (GPU). We demonstrate this process by simultaneously inferring the ice thickness distribution, ice flow parametrization and ice surface of ten of the largest glaciers in Switzerland. As a result, we obtain a high degree of assimilation while guaranteeing an equilibrium between mass-balance and ice flow mechanics. The code runs very efficiently (optimizing one large-size glacier at 100 m takes < 1 min on a laptop) while it is open-source and publicly available.

This paper introduces the Instructed Glacier Model (IGM) – a model that simulates ice dynamics, mass balance and its coupling to predict the evolution of glaciers, icefields or ice sheets. The novelty of IGM is that it models the ice flow by a Convolutional Neural Network, which is trained from data generated with hybrid SIA + SSA or Stokes ice flow models. By doing so, the most computationally demanding model component is substituted by a cheap emulator. Once trained with representative data, we demonstrate that IGM permits to model mountain glaciers up to 1000 × faster than Stokes ones on Central Processing Units (CPU) with fidelity levels above 90% in terms of ice flow solutions leading to nearly identical transient thickness evolution. Switching to the GPU often permits additional significant speed-ups, especially when emulating Stokes dynamics or/and modelling at high spatial resolution. IGM is an open-source Python code which deals with two-dimensional (2-D) gridded input and output data. Together with a companion library of trained ice flow emulators, IGM permits user-friendly, highly efficient and mechanically state-of-the-art glacier and icefields simulations.

Retreat of calving glaciers worldwide has contributed substantially to sea-level rise in recent decades. Mass loss by calving contributes significantly to the uncertainty of sea-level rise projections. At Bowdoin Glacier, Northwest Greenland, most calving occurs by a few large events resulting from kilometre-scale fractures forming parallel to the calving front. High-resolution terrestrial radar interferometry data of such an event reveal that crevasse opening is fastest at low tide and accelerates during the final 36 h before calving. Using the ice flow model Elmer/Ice, we identify the crevasse water level as a key driver of modelled opening rates. Sea water-level variations in the range of local tidal amplitude (1 m) can reproduce observed opening rate fluctuations, provided crevasse water level is at least 4 m above the low-tide sea level. The accelerated opening rates within the final 36 h before calving can be modelled by additional meltwater input into the crevasse, enhanced ice cliff undercutting by submarine melt, ice damage increase due to tidal cyclic fatigue, crevasse deepening or a combination of these processes. Our results highlight the influence of surface meltwater and tides on crevasse opening leading to major calving events at grounded tidewater glaciers such as Bowdoin.

Modelling paleo-glacier networks in mountain ranges on the millennial timescales requires ice flow approximations. Hybrid models calculating ice flow by combining vertical shearing (shallow ice approximation) and longitudinal stretching (shallow shelf approximation) have been applied to model paleo-glacier networks on steep terrain, yet their validity has not yet been assessed quantitatively. Moreover, hybrid models consistently yield higher ice thicknesses than Last Glacial Maximum geomorphological reconstructions in the European Alps. Here, we compare results based on the hybrid Parallel Ice Sheet Model (PISM) and the Stokes model Elmer/Ice on the Rhine Glacier, a catchment of the former European Alpine Icefield. For PISM, we also test two magnitudes of flux limitation in a scheme that reduces shearing velocities. We find that the flux limitation typically used in PISM yields significantly reduced shearing speeds and increases ice thicknesses by up to 500 m, partly explaining previous overestimations. However, reducing the ice flux limitation allows the hybrid model to minimize this mismatch and captures sliding speeds, ice thicknesses, ice extent and basal temperatures in close agreement with those obtained with the Stokes model.

We model the future evolution of the largest glacier of the European Alps – Great Aletsch Glacier, Switzerland – during the 21st century. For that purpose we use a detailed three-dimensional model, which combines full Stokes ice dynamics and surface mass balance forced with the most recent climate projections (CH2018), as well as with climate data of the last decades. As a result, all CH2018 climate scenarios yield a major glacier retreat: Results range from a loss of 60% of today's ice volume by 2100 for a moderate CO2 emission scenario (RCP2.6) being in line with the Paris agreement to an almost complete wastage of the ice for the most extreme emission scenario (RCP8.5). Our model results also provide evidence that half of the mass loss is already committed under the climate conditions of the last decade.

For more than a century Alpine glaciers have been retreating dramatically, and they are expected to shrink even more quickly over the coming decades. This study addresses the future evolution of Grosser Aletschgletscher, Switzerland, the largest glacier in the European Alps. A three-dimensional combined surface mass-balance and glacier dynamics model was applied. The ice flow was described with the full Stokes equations. The glacier surface evolution was obtained by solving a transport equation for the volume of fluid. Daily surface melt and accumulation were calculated on the basis of climate data. The combined model was validated against several types of measurements made throughout the 20th century. For future climate change, scenarios based on regional climate models in the ENSEMBLES project were used. According to the median climatic evolution, Aletschgletscher was expected to lose 90% of its ice volume by the end of 2100. Even when the model was driven using current climate conditions (the past two decades) the glacier tongue experienced a considerable retreat of 6 km, indicating its strong disequilibrium with the present climate. By including a model for the evolution of supraglacial debris and its effect in reducing glacier melt, we show that this factor can significantly slow future glacier retreat.

We propose a novel Eulerian algorithm to compute the changes of a glacier geometry for given mass balances. The surface of a glacier is obtained by solving a transport equation for the volume of fluid (VOF). The surface mass balance is taken into account by adding an interfacial term in the transport equation. An unstructured mesh with standard stabilized finite elements is used to solve the non-linear Stokes problem. The VOF function is computed on a structured grid with a high resolution. The algorithm is stable for Courant numbers larger than unity and conserves mass to a high accuracy. To demonstrate the potential of the algorithm, we apply it to reconstructed Late-glacial states of a small valley glacier, Vadret Muragl, in the Swiss Alps.

The existence of solutions of the non-sliding shallow-ice-sheet equation on a flat horizontal bed with a mass balance linearly depending on altitude is proven for fixed margins. Free-margin solutions for the same mass balance do not exist. Fixed-margin solutions show unbounded shear stress and nonzero mass flux at the margin. Steady-state solutions with realistic margins, vanishing ice flux and vanishing shear stress are found numerically for ice sheets with Weertman-type sliding.

Predictions of marine ice-sheet behaviour require models able to simulate grounding-line migration. We present results of an intercomparison experiment for plan-view marine ice-sheet models. Verification is effected by comparison with approximate analytical solutions for flux across the grounding line using simplified geometrical configurations (no lateral variations, no buttressing effects from lateral drag). Perturbation experiments specifying spatial variation in basal sliding parameters permitted the evolution of curved grounding lines, generating buttressing effects. The experiments showed regions of compression and extensional flow across the grounding line, thereby invalidating the boundary layer theory. Steady-state grounding-line positions were found to be dependent on the level of physical model approximation. Resolving grounding lines requires inclusion of membrane stresses, a sufficiently small grid size (<500 m), or subgrid interpolation of the grounding line. The latter still requires nominal grid sizes of <5 km. For larger grid spacings, appropriate parameterizations for ice flux may be imposed at the grounding line, but the short-time transient behaviour is then incorrect and different from models that do not incorporate grounding-line parameterizations. The numerical error associated with predicting grounding-line motion can be reduced significantly below the errors associated with parameter ignorance and uncertainties in future scenarios.

In this paper we reconstruct the space–time trajectory beneath the surface of Aletschgletscher, Switzerland, of the corpses of three mountaineers that disappeared in March 1926 and reappeared at the glacier surface in June 2012. Our method integrates the time-dependent velocity field of an existing full-Stokes glacier model, starting at the point where the corpses were found at the glacier surface. Our main result is that we were able to localize the immersion location where the brothers presumably died. As a second result, the upstream end point of the computed trajectory emerges very close to the glacier surface in 1926, giving a new and global validation of the glacier model in space and time. Testing the sensitivity of the immersion location obtained with respect to the model and other uncertainties indicates an area of 0.6% of the entire glacier area where the accident could have occurred. Our result suggests that death was not caused by an avalanche or a fall into a crevasse; instead, it is likely that the mountaineers became disoriented in prolonged severe weather conditions and froze to death.

In this study, a modelling approach was used to investigate the cause of the diversion of erratic boulders from Mont Blanc and southern Valais by the Valais Glacier to the Solothurn lobe during the Last Glacial Maximum (LGM). Using the Parallel Ice Sheet Model, we simulated the ice flow field during the LGM, and analyzed the trajectories taken by erratic boulders from areas with characteristic lithologies. The main difficulty in this exercise laid with the large uncertainties affecting the paleo climate forcing required as input for the surface mass-balance model. In order to mimic the prevailing climate conditions during the LGM, we applied different temperature offsets and regional precipitation corrections to present-day climate data, and selected the parametrizations, which yielded the best match between the modelled ice extent and the geomorphologically-based ice-margin reconstruction. After running a range of simulations with varying parameters, our results showed that only one parametrization allowed boulders to be diverted to the Solothurn lobe during the LGM. This precipitation pattern supports the existing theory of preferential southwesterly advection of moisture to the alps during the LGM, but also indicates strongly enhanced precipitation over the Mont Blanc massif and enhanced cooling over the Jura Mountains.

A new hybrid model for the dynamics of glaciers, ice sheets and ice shelves is introduced. In this ‘multilayer’ model the domain of ice consists of a pile of thin layers, which can spread out, contract and slide over each other such that the two most relevant types of stresses are accounted for: membrane and vertical shear. Assuming the horizontal velocity field to be vertically piecewise constant in each layer, the model is derived from local depth integrations of the hydrostatic approximation of the Stokes equations. These integrations give rise to interlayer tractions, which can be redefined at zeroth order in the interlayer surface slope by keeping the vertical shear stress components. Furthermore, if the layers are chosen such that they are aligned with the streamlines, then second-order accurate interlayer tractions can replace zeroth-order ones. The final model consists of a tridiagonal system of two-dimensional nonlinear elliptic equations, the size of this system being equal to the number of layers. When running the model for prognostic flowline ISMIP-HOM benchmark experiments, the multilayer solutions show good agreement with the higher-order solutions if no severe depression occurs in the bedrock. As an alternative to three-dimensional models, the multilayer approach offers to glacier and ice sheet modellers a way of upgrading the commonly used shallow shelf approximation model into a mechanically complete but mathematically two-dimensional model.