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

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.

Free–Surface Flow as a Variational Inequality (<i>evolve_glacier v1.1</i>): Numerical Aspects of a Glaciological Application

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) can not 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. To solve this problem, the 'constrained' non--linear problem solving capabilities of PETSc's SNES 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 no--negativity of ice thickness as well as mass conservation. The simulation framework provides the possibility to partition the computational domain 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 glacierized and ice-free regions. The option to chose different time discretizations, spatial stabilisation schemes and adaptive mesh refinement make it a versatile tool for glaciological applications. We present a set of benchmark tests that highlight 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 discretization methods. We discuss their applicability and suitability for different glaciological applications.

Interaction of a mantle plume and a moving plate: insights from numerical modeling

&lt;p&gt;The impingement of a hot buoyant mantle plume onto the lithosphere can result in either breaking of the lithosphere, which might results in subduction initiation or in under-plating of the plume beneath the lithosphere. Key natural examples of the former and latter are formation of subduction along the southern margin of Caribbean and northwestern South America in the late Cretaceous as well as the hotspot chains of Hawaii, respectively. In previous studies the interaction of a buoyant mantle plume with lithosphere was investigated either for the case of stationary lithosphere or for moving lithosphere but ignoring the effect of magmatic weakening of the lithosphere above the plume head. In this study we aim to investigate the response of a moving lithosphere to the arrival of a stationary mantle plume including the effect of magmatic lithospheric weakening. To do so we use 3d thermo-mechanical models employing the finite difference code I3ELVIS. Our setup consists of an oceanic lithosphere, mantle plume and asthenosphere till depth of 400 km. The moving plate is simulated by imposing a kinematic boundary condition on the lithospheric part of the side boundaries. The mantle plume in our models has a mushroom shape. The experiments differ in the age of the lithosphere, rate of the plate motion and size of the mantle plume. For different combinations of these parameters model results show either (1) breaking of the lithosphere and initiation of subduction above the plume head or (2) asymmetric spreading of the plume material below the lithosphere without large deformation of the lithosphere. We find that the critical radius of the plume that breaks the lithosphere and initiates subduction depends on plume buoyancy and the lithospheric age, but not on the plate speed. In general, the modeling results for the moving plate are similar to the results for a stationary plate, but the shapes of the region of the deformed lithosphere differ.&lt;/p&gt;

Mechanical effect of mélange-induced buttressing on embayment-terminating glacier dynamics

Embayment terminating glaciers interact dynamically with seasonal sea ice and icebergs, a mixture we refer to as mélange. For certain glaciers, mélange prevents calved bergs from rotating away from the front, thus allowing the ice front to advance into the embayment. Here we demonstrate that mélange can, if rigid enough, provide sufficient buttressing to reduce the calving rate, while leaving the ice-front velocity largely unaffected. The net result is additional ice-front advance. Observations indicate a seasonal advance/retreat cycle has occurred at Jakobshavn Isbræ since the 1950s. We model an idealized Jakobshavn Isbræ-like scenario and find that mélange may be responsible for a seasonal ice-front advance of up to 0.6 km. These results come from a model that incorporates mélange into the interior of the domain, includes relevant stresses, and models drag via a kinematic boundary condition. A weakening or loss of mélange due to increasing temperatures would lead to further mass loss from glaciers such as Jakobshavn Isbræ.