Inequality-constrained free-surface evolution in a full Stokes ice flow model (<i>evolve_glacier v1.1</i>)

Abstract

Abstract. Like any gravitationally driven flow that is not constrained at the upper surface, glaciers and ice sheets feature a free surface, which becomes a free-boundary problem within simulations. A kinematic boundary condition is often used to describe the evolution of this free surface. However, in the case of glaciers and ice sheets, the naturally occurring constraint that the ice surface elevation (S) cannot fall below the bed topography (B) (S-B≥0), in combination with a non-zero mass balance rate complicates the matter substantially. We present an open-source numerical simulation framework to simulate the free-surface evolution of glaciers that directly incorporates this natural constraint. It is based on the finite-element software package FEniCS solving the Stokes equations for ice flow and a suitable transport equation, i.e. “kinematic boundary condition”, for the free-surface evolution. The evolution of the free surface is treated as a variational inequality, constrained by the bedrock underlying the glacier or the topography of the surrounding ground. This problem is solved using a “reduced space” method, where a Newton line search is performed on a subset of the problem (Benson and Munson, 2006). Therefore, the “constrained” non-linear problem-solving capabilities of PETSc's (Portable, Extensible Toolkit for Scientific Computation, Balay et al., 2019) SNES (Scalable Non-linear Equations Solver) interface are used. As the constraint is considered in the solving process, this approach does not require any ad hoc post-processing steps to enforce non-negativity of ice thickness and corresponding mass conservation. The simulation framework provides the possibility to divide the computational domain into different subdomains so that individual forms of the relevant equations can be solved for different subdomains all at once. In the presented setup, this is used to distinguish between glacierised and ice-free regions. The option to chose different time discretisations, spatial stabilisation schemes and adaptive mesh refinement make it a versatile tool for glaciological applications. We present a set of benchmark tests that highlight that the simulation framework is able to reproduce the free-surface evolution of complex geometries under different conditions for which it is mass-conserving and numerically stable. Real-world glacier examples demonstrate high-resolution change in glacier geometry due to fully resolved 3D velocities and spatially variable mass balance rate, whereby realistic glacier recession and advance states can be simulated. Additionally, we provide a thorough analysis of different spatial stabilisation techniques as well as time discretisation methods. We discuss their applicability and suitability for different glaciological applications.