Shallow Ice Approximation Research Articles

Surface Motion and Topographic Effects on Ice Thickness Inversion for High Mountain Asia Glaciers: A Comparison Study from Three Numerical Models

Several studies have evaluated the accuracy of ice thickness model inversions, but the effects of surface velocity and topographic factors have not been fully studied. Because of the important potential of models for estimating regional and global ice thickness, the influence of the characteristics of the glaciers themselves on the models is unknown. Therefore, we estimated the thickness of 10 glaciers in HMA based on remote sensing data using the GlabTop2, GV22, and VOLTA models, and we analyzed the factors affecting accuracy. The first two models are representative of the widely used shallow ice approximation and laminar flow theory, and the latter is a new model based on shallow ice approximation. The results show that the RMSEs for GlabTop2, VOLTA, and GV22 are 36.6 m, 56.8 m, and 63.2 m, respectively. The accuracy of the three models differed obviously under the influence of different slopes, aspects, and surface velocities. In contrast to glaciers with higher surface velocities, GV22 is more suitable for glaciers with smaller surface velocities, and the accuracy of the model is significantly reduced when surface velocities are greater than 50 m/yr. Aspect has an effect on the accuracy of the GV22 ice thickness inversion but not on the other models. Ice thicknesses estimated using the three models is more reliable at surface slopes of 16–24°. These three models do not apply to gentle glaciers such as icecaps. Our study provides an important reference for choosing a suitable model to invert the glacier thickness in HMA.

Universal differential equations for glacier ice flow modelling

Abstract. Geoscientific models are facing increasing challenges to exploit growing datasets coming from remote sensing. Universal differential equations (UDEs), aided by differentiable programming, provide a new scientific modelling paradigm enabling both complex functional inversions to potentially discover new physical laws and data assimilation from heterogeneous and sparse observations. We demonstrate an application of UDEs as a proof of concept to learn the creep component of ice flow, i.e. a nonlinear diffusivity differential equation, of a glacier evolution model. By combining a mechanistic model based on a two-dimensional shallow-ice approximation partial differential equation with an embedded neural network, i.e. a UDE, we can learn parts of an equation as nonlinear functions that then can be translated into mathematical expressions. We implemented this modelling framework as ODINN.jl, a package in the Julia programming language, providing high performance, source-to-source automatic differentiation (AD) and seamless integration with tools and global datasets from the Open Global Glacier Model in Python. We demonstrate this concept for 17 different glaciers around the world, for which we successfully recover a prescribed artificial law describing ice creep variability by solving ∼ 500 000 ordinary differential equations in parallel. Furthermore, we investigate which are the best tools in the scientific machine learning ecosystem in Julia to differentiate and optimize large nonlinear diffusivity UDEs. This study represents a proof of concept for a new modelling framework aiming at discovering empirical laws for large-scale glacier processes, such as the variability in ice creep and basal sliding for ice flow, and new hybrid surface mass balance models.

Modeling large‐scale landform evolution with a stream power law for glacial erosion (OpenLEM v37): benchmarking experiments against a more process-based description of ice flow (iSOSIA v3.4.3)

Abstract. Following the tradition of modeling fluvial landscape evolution, a novel approach describing glacial erosion based on an empirical stream power law was proposed. This approach differs substantially from well-established process-based models applied to describe glacial erosion in mountain landscapes. Outstanding computational performance but a number of potential limitations compared to process-based models requires extensive testing to evaluate the applicability of this novel approach. In this study, we test the validity of the glacial stream power law and its implementation into a 2-D landform evolution model (OpenLEM) by benchmarking it against a state of the art surface process model based on the integrated second-order shallow-ice approximation (iSOSIA). Despite completely different approaches, OpenLEM and iSOSIA predict similar ice flow patterns and erosion rates for a wide range of climatic conditions without re-adjusting a set of calibrated parameters. This parameter set is valid for full glacial conditions where the entire precipitation is converted to ice but also for an altitude-dependent glacier mass balance as characteristic for most glaciated mountain ranges on Earth. In both models characteristic glacial features, such as overdeepenings, hanging valleys and steps at confluences emerge roughly at the same locations, resulting in a consistent altitude-dependent adjustment of channel slope and relief. Compared to iSOSIA, however, distinctly higher erosion rates occur in OpenLEM at valley flanks during the initial phase of the fluvial to glacial transition. This is mainly due to the simplified description of glacier width and ice surface in OpenLEM. In this respect, we found that the glacial stream power approach cannot replace process-based models such as iSOSIA but is complementary to them by addressing research questions that could not previously be answered due to a lack of computational efficiency. The implementation of the glacial stream power law is primarily suitable for large-scale simulations investigating the evolution of mountain topography in the interplay of tectonics and climate. As coupling glacial and fluvial erosion with sediment transport shows nearly the same computational efficiency as its purely fluvial counterpart, mountain-range-scale simulations at high spatial resolution are not exclusively restricted to the fluvial domain anymore, and a series of exciting research questions can be addressed by this novel approach.

Performance analysis of high-resolution ice-sheet simulations

AbstractNumerical glacier and ice-sheet models compute evolving ice geometry and velocity fields using various stress-balance approximations and boundary conditions. At high spatial resolution, with horizontal mesh/grid resolutions of a few kilometers or smaller, these models usually require time steps shorter than climate-coupling time scales because they update ice thickness after each velocity solution. High-resolution performance is degraded by the stability restrictions of such explicit time-stepping. This short note, which considers the shallow ice approximation and Stokes models as stress-balance end members, clarifies the scaling of numerical model performance by quantifying simulation cost per model year in terms of mesh resolution and the number of degrees of freedom. The performance of current-generation explicit time-stepping models is assessed, and then compared to the prospective performance of implicit schemes. The main results highlight the key roles played by the algorithmic scaling of stress-balance solvers and coupled, implicit-step solvers.

Hysteretic evolution of ice rises and ice rumples in response to variations in sea level

Abstract. Ice rises and ice rumples are locally grounded features found in coastal Antarctica and are surrounded by otherwise freely floating ice shelves. An ice rise has an independent flow regime, whereas the flow regime of an ice rumple conforms to that of the ice shelf and merely slows the flow of ice. In both cases, local highs in the bathymetry are in contact with the ice shelf from below, thereby regulating the large-scale ice flow, with implications for the upstream continental grounding line position. This buttressing effect, paired with the suitability of ice rises as a climate archive, necessitates a better understanding of the transition between ice rise and ice rumple, their evolution in response to a change in sea level, and their dynamic interaction with the surrounding ice shelf. We investigate this behaviour using a three-dimensional full Stokes ice flow model with idealised ice rises and ice rumples. The simulations span end-member basal friction scenarios of almost stagnant and fully sliding ice at the ice–bed interface. We analyse the coupling with the surrounding ice shelf by comparing the deviations between the non-local full Stokes surface velocities and the local shallow ice approximation (SIA). Deviations are generally high at the ice divides and small on the lee sides. On the stoss side, where ice rise and ice shelf have opposing flow directions, deviations can be significant. Differences are negligible in the absence of basal sliding where the corresponding steady-state ice rise is larger and develops a fully independent flow regime that is well described by SIA. When sea level is increased, and a transition from ice rise to ice rumple is approached, the divide migration is more abrupt the higher the basal friction. In each scenario, the transition occurs after the stoss-side grounding line has moved over the bed high and is positioned on a retrograde slope. We identify a hysteretic response of ice rises and ice rumples to changes in sea level, with grounded area being larger in a sea-level-increase scenario than in a sea-level-decrease scenario. This hysteresis shows not only irreversibility following an equal increase and subsequent decrease in sea level but also that the perturbation history is important when the ice rise or ice rumple geometry is not known. The initial grounded area needs to be carefully considered, as this will determine the formation of either an ice rise or an ice rumple, thereby causing different buttressing effects.

Inferring rheological properties and topographical features from free surface flow data

The accurate modelling of geophysical flows often requires information that is difficult to measure and therefore poorly quantified. Such information may relate to the fluid properties or an unknown boundary condition, for example. The premise of this talk is that when the flow is bounded by a free surface, the deformation of this free surface contains useful information, which can be used to infer such unknown quantities. The increasing availability of free surface data through remote sensing using drones and satellites provides the impetus to use mathematical methods and numerical tools to interpret the signature embedded in the free surface deformation.In the first part of the talk, we will explore the problem of recovering simultaneously the ice thickness and basal slip of an ice flow governed by the shallow ice approximation. In the second part, using similar techniques, we will show how we can recover rheological parameters for thin-film lava flow.

A note on the snout

AbstractThe shallow ice approximation for glaciers and ice sheets is a degenerate model in which the ice surface slope at the margin may be infinite. This result is due to the neglect of the otherwise small longitudinal stress terms. Here we derive a corrected approximation for the basal shear stress, and show that the resulting model provides an explanation for the observed finite slope margins.

GREB-ISM v1.0: A coupled ice sheet model for the Globally Resolved Energy Balance model for global simulations on timescales of 100 kyr

Abstract. We introduce a newly developed global ice sheet model coupled to the Globally Resolved Energy Balance (GREB) climate model for the simulation of global ice sheet evolution on timescales of 100 kyr or longer (GREB-ISM v1.0). Ice sheets and ice shelves are simulated on a global grid, fully interacting with the climate simulation of surface temperature, precipitation, albedo, land–sea mask, topography and sea level. Thus, it is a fully coupled atmosphere, ocean, land and ice sheet model. We test the model in ice sheet stand-alone and fully coupled simulations. The ice sheet model dynamics behave similarly to other hybrid SIA (shallow ice approximation) and SSA (shallow shelf approximation) models, but the West Antarctic Ice Sheet accumulates too much ice using present-day boundary conditions. The coupled model simulations produce global equilibrium ice sheet volumes and calving rates like those observed for present-day boundary conditions. We designed a series of idealized experiments driven by oscillating solar radiation forcing on periods of 20, 50 and 100 kyr in the Northern Hemisphere. These simulations show clear interactions between the climate system and ice sheets, resulting in slow buildup and fast decay of ice-covered areas and global ice volume. The results also show that Northern Hemisphere ice sheets respond more strongly to timescales longer than 100 kyr. The coupling to the atmosphere and sea level leads to climate interactions between the Northern and Southern Hemispheres. The model can run global simulations of 100 kyr d−1 on a desktop computer, allowing the simulation of the whole Quaternary period (2.6 Myr) within 1 month.

Brief communication: Estimating the ice thickness of the Müller Ice Cap to support selection of a drill site

Abstract. The Müller Ice Cap will soon set the scene for a new drilling project. Therefore, ice thickness estimates are necessary for planning, since thickness measurements of the ice cap are sparse. Here, three models are presented and compared: (i) a simple Semi-Empirical Ice Thickness Model (SEITMo) based on an inversion of the shallow-ice approximation by the use of a single radar line in combination with the glacier outline, surface slope, and elevation; (ii) an iterative inverse method using the Parallel Ice Sheet Model (PISM), and (iii) a velocity-based inversion of the shallow-ice approximation. The velocity-based inversion underestimates the ice thickness at the ice cap top, making the model less useful to aid in drill site selection, whereas PISM and the SEITMo mostly agree about a good drill site candidate. However, the new SEITMo is insensitive to mass balance, computationally fast, and provides as good fits as PISM.

An augmented lagrangian algorithm for recovery of ice thickness in unidirectional flow using the shallow ice approximation

A key parameter in ice flow modelling is basal slipping at the ice-bed interface as it can have a large effect on the resultant ice thickness. Unfortunately, its contribution to surface observations can be hard to distinguish from that of bed undulations. Therefore, inferring the ice thickness from surface measurements is an interesting and non-trivial inverse problem. This paper presents a method for recovering dually the ice thickness, in turn the bed topography, and the basal slip using only surface elevation and speed measurements. The unidirectional shallow ice approximation is first implemented to model steady state ice flow for given bedrock and basal slip profiles. This surface is then taken as synthetic observed data. To invert the data, an optimisation approach is taken, employing a basic Uzawa type algorithm to minimise the augmented Lagrangian objective functional. This optimisation recovers the diffusion coefficient which gives the best fit to observations. Combining this recovered diffusion with observed surface velocity, a simple Newton’s method is used to recover both the ice thickness and basal slip. The method was successful in each test case and this implies that it should be possible to recover both of these parameters in two-dimensional cases also.

Physically-constrained data-driven inversions to infer the bed topography beneath glaciers flows. Application to East Antarctica

A method for infering bed topography beneath glaciers from surface measurements (elevation from altimetry and velocity from InSAR) and sparse thickness measurements is developed and evaluated. The method is based on an original non-isothermal Reduced Uncertainty (RU) version of the Shallow Ice Approximation (SIA) equation that natively incorporates the surface measurements. The flow model has a single dimensionless multi-physics parameter γ. This parameter takes into account the basal slipperiness and the variable vertical rate factor profiles, thus the vertical thermal variations. The inversions are based on three steps involving: an Artificial Neural Network (ANN) and two Variational Data Assimilation (VDA) processes. The ANN-based stage aims at estimating the multi-physics number γ from the thickness measurements; the resulting estimator is remarkably robust. The full inversion method is valid for half-sheared flows (presenting a moderate basal slipperiness): it can be applied to inland ice-sheets areas. Also these estimates connect continuously with estimates from mass conservation only, i.e. with areas of sliding flows. Numerical results are presented for areas of the East Antarctica Ice Sheet where bed elevation can be very uncertain (Bedmap2 values). Estimates are valid for wavelengths longer than \(\sim 10 \bar h\) (due to the long wave assumption, shallow flow model) with resolution at \(\sim \bar h\) (\(\bar h\) a characteristic thickness value).

A Higher-Order Singular Value Decomposition Tensor Emulator for Spatiotemporal Simulators

We introduce methodology to construct an emulator for environmental and ecological spatiotemporal processes that uses the higher-order singular value decomposition (HOSVD) as an extension of singular value decomposition (SVD) approaches to emulation. Some important advantages of the method are that it allows for the use of a combination of supervised learning methods (e.g., random forests and Gaussian process regression) and also allows for the prediction of process values at spatial locations and time points that were not used in the training sample. The method is demonstrated with two applications: The first is a periodic solution to a shallow ice approximation partial differential equation from glaciology, and second is an agent-based model of collective animal movement. In both cases, we demonstrate the value of combining different machine learning models for accurate emulation. In addition, in the agent-based model case we demonstrate the ability of the tensor emulator to successfully capture individual behavior in space and time. We demonstrate via a real data example the ability to perform Bayesian inference in order to learn parameters governing collective animal behavior.

Bayesian Inference of Ice Softness and Basal Sliding Parameters at Langjökull

We develop Bayesian statistical models that are designed for the inference of ice softness and basal sliding parameters, important glaciological quantities. These models are applied to Langjökull, the second largest temperate ice cap in Iceland at about 900 squared kilometers in area. The models make use of a relationship between physical parameters and ice velocity as stipulated by a shallow ice approximation that is generally applicable to Langjökull. The posterior distribution for ice softness concentrates around 18.2 × 10−25s−1Pa−3; moreover, spatially varying basal sliding parameters are inferred allowing for the decomposition of velocity into a deformation component and a sliding component, with spatial variation consistent with previous studies. Bayesian computation is conducted with a Gibbs sampling approach. The paper serves as an example of statistical inference for ice softness and basal sliding parameters at temperate, shallow glaciers using surface velocity data.

Simulation, modeling and authoring of glaciers

Glaciers are some of the most visually arresting and scenic elements of cold regions and high mountain landscapes. Although snow-covered terrains have previously received attention in computer graphics, simulating the temporal evolution of glaciers as well as modeling their wide range of features has never been addressed. In this paper, we combine a Shallow Ice Approximation simulation with a procedural amplification process to author high-resolution realistic glaciers. Our multiresolution method allows the interactive simulation of the formation and the evolution of glaciers over hundreds of years. The user can easily modify the environment variables, such as the average temperature or precipitation rate, to control the glacier growth, or directly use brushes to sculpt the ice or bedrock with interactive feedback. Mesoscale and smallscale landforms that are not captured by the glacier simulation, such as crevasses, moraines, seracs, ogives, or icefalls are synthesized using procedural rules inspired by observations in glaciology and according to the physical parameters derived from the simulation. Our method lends itself to seamless integration into production pipelines to decorate reliefs with glaciers and realistic ice features.

Bedrock reconstruction from free surface data for unidirectional glacier flow with basal slip

Glacier ice flow is shaped and defined by several properties, including the bedrock elevation profile and the basal slip distribution. The effect of these two basal properties can present in similar ways in the surface. For bedrock recovery this makes distinguishing between them an interesting and complex problem. The results of this paper show that in some synthetic test cases it is indeed possible to distinguish and recover both bedrock elevation and basal slip given free surface elevation and free surface velocity. The unidirectional shallow ice approximation is used to compute steady state surface data for a number of synthetic cases with different bedrock profiles and basal slip distributions. A simple inversion method based on Newton's method is applied to the known surface data to return the bedrock profile and basal slip distribution. In each synthetic test case, the inversion was successful in recovering both the bedrock elevation profile and the basal slip distribution variables. These results imply that there is a unique bedrock profile and basal slip which give rise to a unique combination of free surface velocity and free surface elevation.

Lagrangian formulation, a general relativity analogue, and a symmetry of the Vialov equation of glaciology

Using a suitable rescaling of the independent variable, a Lagrangian is found for the nonlinear Vialov equation ruling the longitudinal profiles of glaciers and ice caps in the shallow ice approximation. This leads to a formal analogy between the (rescaled) Vialov equation and the Friedmann equation of relativistic cosmology, which is explored. This context provides a new symmetry of the (rescaled) Vialov equation and gives, at least formally, all its solutions using a generating function, which is the Nye profile for the degenerate case of perfectly plastic ice.

Possible biases in scaling-based estimates of glacier change: a case study in the Himalaya

Abstract. Approximate glacier models are routinely used to compute the future evolution of mountain glaciers under any given climate-change scenario. A majority of these models are based on statistical scaling relations between glacier volume, area, and/or length. In this paper, long-term predictions from scaling-based models are compared with those from a two-dimensional shallow-ice approximation (SIA) model. We derive expressions for climate sensitivity and response time of glaciers assuming a time-independent volume–area scaling. These expressions are validated using a scaling-model simulation of the response of 703 synthetic glaciers from the central Himalaya to a step change in climate. The same experiment repeated with the SIA model yields about 2 times larger climate sensitivity and response time than those predicted by the scaling model. In addition, the SIA model obtains area response time that is about 1.5 times larger than the corresponding volume response time, whereas scaling models implicitly assume the two response times to be equal to each other. These results indicate the possibility of a low bias in the scaling model estimates of the long-term loss of glacier area and volume. The SIA model outputs are used to obtain parameterisations, climate sensitivity, and response time of glaciers as functions of ablation rate near the terminus, mass-balance gradient, and mean thickness. Using a linear-response model based on these parameterisations, we find that the linear-response model outperforms the scaling model in reproducing the glacier response simulated by the SIA model. This linear-response model may be useful for predicting the evolution of mountain glaciers on a global scale.

A glacial buzzsaw effect generated by efficient erosion of temperate glaciers in a steady state model

Studies on relief and hypsometry in glaciated mountain ranges have recognized the equilibrium line altitude of glaciers as a major control on mountain height. This has led to the belief that, in contrast to fluvial topography, glacial topography is limited by climate, with tectonics playing a secondary role. This concept is known as the glacial buzzsaw. However, the understanding of controls on glacial relief has remained mainly qualitative, in part because a reference scenario for glacial landscape evolution — a glacial topographic steady state linked to the base level by a fluvial topographic steady-state downstream of the glacier terminus - and its dependencies have not yet been defined. Here we define such a reference and compare steady state longitudinal profiles in a coupled system of glacial and fluvial erosion that involves both tectonics and climate. Our model is based on the stream power law for fluvial erosion, the shallow ice approximation in combination with a glacial erosion rule, and an empirically determined drainage area-flow length relationship for both rivers and glaciers. Further, we introduce a new approach to incorporate dendritic glacier network structures into the one-dimensional model. Modeling of coupled glacial-fluvial steady-state profiles with empirical glacial and fluvial erosion parameters shows that the difference between glacial and fluvial relief depends on climate, tectonics and the applied erosion laws. Our results imply that glacial erosion can typically balance tectonic uplift rates at lower relief and topographic slopes than fluvial erosion. This steady state equivalent of the glacial buzzsaw effect suggests that glaciers may indeed be more efficient erosional agents than rivers increasing erosion in a cooling climate.

Modelling a paleo valley glacier network using a hybrid model: an assessment with a Stokes ice flow model

AbstractModelling 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.

The Open Global Glacier Model (OGGM) v1.1

Abstract. Despite their importance for sea-level rise, seasonal water availability, and as a source of geohazards, mountain glaciers are one of the few remaining subsystems of the global climate system for which no globally applicable, open source, community-driven model exists. Here we present the Open Global Glacier Model (OGGM), developed to provide a modular and open-source numerical model framework for simulating past and future change of any glacier in the world. The modeling chain comprises data downloading tools (glacier outlines, topography, climate, validation data), a preprocessing module, a mass-balance model, a distributed ice thickness estimation model, and an ice-flow model. The monthly mass balance is obtained from gridded climate data and a temperature index melt model. To our knowledge, OGGM is the first global model to explicitly simulate glacier dynamics: the model relies on the shallow-ice approximation to compute the depth-integrated flux of ice along multiple connected flow lines. In this paper, we describe and illustrate each processing step by applying the model to a selection of glaciers before running global simulations under idealized climate forcings. Even without an in-depth calibration, the model shows very realistic behavior. We are able to reproduce earlier estimates of global glacier volume by varying the ice dynamical parameters within a range of plausible values. At the same time, the increased complexity of OGGM compared to other prevalent global glacier models comes at a reasonable computational cost: several dozen glaciers can be simulated on a personal computer, whereas global simulations realized in a supercomputing environment take up to a few hours per century. Thanks to the modular framework, modules of various complexity can be added to the code base, which allows for new kinds of model intercomparison studies in a controlled environment. Future developments will add new physical processes to the model as well as automated calibration tools. Extensions or alternative parameterizations can be easily added by the community thanks to comprehensive documentation. OGGM spans a wide range of applications, from ice–climate interaction studies at millennial timescales to estimates of the contribution of glaciers to past and future sea-level change. It has the potential to become a self-sustained community-driven model for global and regional glacier evolution.