A comparison of the stability and performance of depth-integrated ice-dynamics solvers

In the last decade, the number of ice-sheet models has increased substantially, in line with the growth of the glaciological community. These models use solvers based on different approximations of ice dynamics. In particular, several depth-integrated dynamics solvers have emerged as fast solvers capable of resolving the relevant physics of ice sheets at the continental scale. However, the numerical stability of these schemes has not been studied systematically to evaluate their effectiveness in practice. Here we focus on three such solvers, the so-called Hybrid, L1L2-SIA and DIVA solvers, as well as the well-known SIA and SSA solvers as boundary cases. We investigate the numerical stability of these solvers as a function of grid resolution and the state of the ice sheet for an explicit time discretization scheme of the mass conservation step. Under simplified conditions with constant viscosity, the maximum stable time step of the Hybrid solver, like the SIA solver, has a quadratic dependence on grid resolution. In contrast, the DIVA solver has a maximum time step that is independent of resolution as the grid becomes increasingly refined, like the SSA solver. A simple 1D implementation of the L1L2-SIA solver indicates that it should behave similarly, but in practice, the complexity of its implementation appears to restrict its stability. In realistic simulations of the Greenland Ice Sheet with a nonlinear rheology, the DIVA and SSA solvers maintain superior numerical stability, while the SIA, Hybrid and L1L2-SIA solvers show markedly poorer performance. At a grid resolution of Delta x = 4 km, the DIVA solver runs approximately 20 times faster than the Hybrid and L1L2-SIA solvers as a result of a larger stable time step. Our analysis shows that as resolution increases, the ice-dynamics solver can act as a bottleneck to model performance. The DIVA solver emerges as a clear outlier in terms of both model performance and its representation of the ice-flow physics itself

The LifeV library: engineering mathematics beyond the proof of concept

LifeV is a library for the finite element (FE) solution of partial differential equations in one, two, and three dimensions. It is written in C++ and designed to run on diverse parallel architectures, including cloud and high performance computing facilities. In spite of its academic research nature, meaning a library for the development and testing of new methods, one distinguishing feature of LifeV is its use on real world problems and it is intended to provide a tool for many engineering applications. It has been actually used in computational hemodynamics, including cardiac mechanics and fluid-structure interaction problems, in porous media, ice sheets dynamics for both forward and inverse problems. In this paper we give a short overview of the features of LifeV and its coding paradigms on simple problems. The main focus is on the parallel environment which is mainly driven by domain decomposition methods and based on external libraries such as MPI, the Trilinos project, HDF5 and ParMetis. Dedicated to the memory of Fausto Saleri.Comment: Review of the LifeV Finite Element librar

Parameterization for subgrid-scale motion of ice-shelf calving fronts

In order to explore the response of the Greenland ice sheet (GIS) to climate change on long (centennial to multi-millennial) time scales, a regional energy-moisture balance model has been developed. This model simulates seasonal variations of temperature and precipitation over Greenland and explicitly accounts for elevation and albedo feedbacks. From these fields, the annual mean surface temperature and surface mass balance can be determined and used to force an ice sheet model. The melt component of the surface mass balance is computed here using both a positive degree day approach and a more physically-based alternative that includes insolation and albedo explicitly. As a validation of the climate model, we first simulated temperature and precipitation over Greenland for the prescribed, present-day topography. Our simulated climatology compares well to observations and does not differ significantly from that of a simple parameterization used in many previous simulations. Furthermore, the calculated surface mass balance using both melt schemes falls within the range of recent regional climate model results. For a prescribed, ice-free state, the differences in simulated climatology between the regional energy-moisture balance model and the simple parameterization become significant, with our model showing much stronger summer warming. When coupled to a three-dimensional ice sheet model and initialized with present-day conditions, the two melt schemes both allow realistic simulations of the present-day GIS

Parameterization for subgrid-scale motion of ice-shelf calving fronts

A parameterization for the motion of ice-shelf fronts on a Cartesian grid in finite-difference land-ice models is presented. The scheme prevents artificial thinning of the ice shelf at its edge, which occurs due to the finite resolution of the model. The intuitive numerical implementation diminishes numerical dispersion at the ice front and enables the application of physical boundary conditions to improve the calculation of stress and velocity fields throughout the ice-sheet-shelf system. Numerical properties of this subgrid modification are assessed in the Potsdam Parallel Ice Sheet Model (PISM-PIK) for different geometries in one and two horizontal dimensions and are verified against an analytical solution in a flow-line setup

Variational Methods in Ice Sheet Modelling

A complete simulation of flowing ice requires knowledge of both the fundamental physical principles that govern the stress and energy balances and a framework for assimilating data into a model to help estimate unknown parameters. Modelling ice is complex, due to the large spatial extent of ice sheets, the multiple scales at which relevant physics operate, and the coupling between heat, stress, and ice rheology. As such, it is usually necessary to make approximations to the equations governing ice flow. At the same time, it is important to have an understanding of the specific assumptions that lead to these approximations. We develop a variational principle for Stokes flow, and neglect certain components in order to obtain the variational principle for the first-order approximation for ice flow. This result is fundamentally the result of assuming bed slopes to be much less than surface slopes, and that vertical resistive stresses are negligible. From a practical standpoint, using automatic differentiation tools on this functional yields a compact model of ice flow that automatically incorporates correct boundary condition. This model is compared to well known benchmark tests. We also present an improved model of ice thermodynamics that operates on enthalpy rather than temperature, avoiding many of the difficulties associated with phase change. We derive a method for inverting the Blatter-Pattyn ice sheet model in order to solve for the rate of basal sliding. This method uses the adjoint equations of the forward model to obtain the gradient of an error functional, and this is minimized using a quasi-Newton method. These methods are applied to an instrumented streamline of the Greenland ice sheet. We perform numerical experiments on this geometry in order to assess the sensitivity of thermal conditions at the ice sheet bed to perturbations in unknown parameters. The basal thermal regime is sensitive to changes in geothermal heat flux, with the location of the transition zone between cold and temperate ice being linear sensitive to changes in it. The temperature field of the ice sheet is insensitive to downstream changes in sliding speed due to the short length scales over which longitudinal coupling acts

Coupled ice shelf-ocean modeling and complex grounding line retreat from a seabed ridge

Recent observations and modeling work have shown a complex mechanical coupling between Antarctica's floating ice shelves and the adjacent grounded ice sheet. A prime example is Pine Island Glacier, West Antarctica, which has a strong negative mass balance caused by a recent increase in ocean-induced melting of its ice shelf. The mass loss coincides with the retreat of the grounding line from a seabed ridge, on which it was at least partly grounded until the 1970s. At present, it is unclear what has caused the onset of this retreat and how feedback mechanisms between the ocean and ice shelf geometry have influenced the ice dynamics. To address these questions, we present the first results from an offline coupling between a state-of-the-art shallow-ice flow model with grounding line resolving capabilities and a three-dimensional ocean general circulation model with a static implementation of the ice shelf. A series of idealized experiments simulate the retreat from a seabed ridge in response to changes in the ocean forcing, and we show that the retreat becomes irreversible after 20 years of warm ocean conditions. A comparison to experiments with a simple depth-dependent melt rate parameterization demonstrates that such parameterizations are unable to capture the details of the retreat process, and they overestimate mass loss by more than 40% over a 50 year timescal

Data Assimilation and Prognostic Whole Ice Sheet Modelling with the Variationally Derived, Higher Order, Open Source, and Fully Parallel Ice Sheet Model VarGlaS

We introduce a novel, higher order, finite element ice sheet model called VarGlaS (Variational Glacier Simulator), which is built on the finite element framework FEniCS. Contrary to standard procedure in ice sheet modelling, Var- GlaS formulates ice sheet motion as the minimization of an energy functional, conferring advantages such as a consistent platform for making numerical approximations, a coherent relationship between motion and heat generation, and implicit boundary treatment. VarGlaS also solves the equations of enthalpy rather than temperature, avoiding the solution of a contact problem. Rather than include a lengthy model spinup procedure, VarGlaS possesses an automated framework for model inversion. These capabilities are brought to bear on several benchmark problems in ice sheet modelling, as well as a 500 yr simulation of the Greenland ice sheet at high resolution. VarGlaS performs well in benchmarking experiments and, given a constant climate and a 100 yr relaxation period, predicts a mass evolution of the Greenland ice sheet that matches present-day observations of mass loss. VarGlaS predicts a thinning in the interior and thickening of the margins of the ice sheet