Multilayer shallow shelf approximation: Minimisation formulation, finite element solvers and applications

Abstract

In this paper, a multilayer generalisation of the Shallow Shelf Approximation (SSA) is considered. In this recent hybrid ice flow model, the ice thickness is divided into thin layers, which can spread out, contract and slide over each other in such a way that the velocity profile is layer-wise constant. Like the SSA (1-layer model), the multilayer model can be reformulated as a minimisation problem. However, unlike the SSA, the functional to be minimised involves a new penalisation term for the interlayer jumps of the velocity, which represents the vertical shear stresses induced by interlayer sliding. Taking advantage of this reformulation, numerical solvers developed for the SSA can be naturally extended layer-wise or column-wise. Numerical results show that the column-wise extension of a Newton multigrid solver proves to be robust in the sense that its convergence is barely influenced by the number of layers and the type of ice flow. In addition, the multilayer formulation appears to be naturally better conditioned than the one of the first-order approximation to face the anisotropic conditions of the sliding-dominant ice flow of ISMIP-HOM experiments.

Introduction

Glaciologists need ice flow models which can be run at very large scales (in space and time) and treat the mechanics adequately while being computationally tractable. Examples of applications are in marine ice sheet modelling in order to better evaluate sea level rise in a climate change regime [34], or in paleoglaciology in order to reconstruct the extents of glaciers and ice sheets over glacial cycles [25]. Despite some recent progress achieved in parallelising solvers [2] or domain decomposition [31], 3D models remain too computationally demanding to be used for that purpose. In addition, remeshing procedures in 3D are complex to apply, since meshes must conserve a certain quality in order to preserve the performances of the solver [16]. For this reason, simplified zero-order models of reduced complexity like the Shallow Ice Approximation (SIA) [12] or the Shallow Shelf Approximation (SSA) [18], [17], [35] are still popular in the community of glacier and ice sheet modellers for running large-scale simulations. Based on the assumption of small aspect ratios of the ice geometries, the SIA, which is a mathematical 1D (vertical) model, accounts only for vertical shear stresses, while the SSA, which is mathematically 2D (horizontal), accounts only for longitudinal stresses. If it is justified to use either the SSA or the SIA in some localised parts of the ice domain, it is often necessary to combine the two when modelling the ice flows of an entire glacier or ice sheet. For example, the vertical shear components of the stress tensor are significant where ice is frozen to the ground, while the longitudinal components are dominant in the fast-sliding parts, like the floating areas, such that using the SIA on the grounded part and the SSA on the floating area might look acceptable. Unfortunately, such an approach is not suitable in the vicinity of the Grounding Line (GL) which delimits the grounded and the floating areas, as all components must be accounted for [29]. This has driven the construction of “hybrid” models, which account for both kind of stresses, while being mathematically 2D. The simplest hybrid model consists of the linear combination SIA+SSA, which is arrived at by adding together the velocities of each model [3]. Unfortunately, this model does not include the simultaneous coupling between the vertical shear and the longitudinal stresses. As a result, the model cannot capture the 3D ice flows that occur in deep and narrow valleys or in the vicinity of GLs [23], [14]. In contrast, the L1L2 [11] or some variants like the ones proposed in [24], [30] or in [8] include the vertical shear stress in the computation of the effective viscosity of the SSA. All these hybrid models have in common is that they solve a single non-linear elliptic 2D problem, and that the velocity profile is reconstructed a posteriori via an implicit relation [30], [36], [4].

In this paper, a recently introduced hybrid multilayer model generalising the SSA is considered [13]. This approach consists of seeing the ice thickness as a pile of thin layers which can spread out, contract and slide over each other. Assuming a vertically piecewise-constant velocity profile in each layer, the model derives from local depth-integrations of the hydrostatic approximation [1], [20]. The crucial step when deriving the model consists of redefining the interlayer tractions by keeping only the shear components. The final multilayer model consists of a tridiagonal system of 2D non-linear elliptic equations, whose number corresponds to the number of layers. By construction, this multilayer model naturally generalises the SSA, which corresponds to the 1-layer case of the model. As a consequence, it is called “Multilayer Shallow Shelf Approximation” in this paper and is abbreviated as MSSA. Unlike the SSA, the MSSA is hybrid since it combines the longitudinal and the vertical shear stresses. Like the SSA, the MSSA model can be reformulated as a minimisation problem for a certain functional. Interestingly, with such a reformulation, the new term corresponding to vertical shear stresses can been interpreted as a penalisation term for the interlayer jumps of the velocity components. Finally, in contrast to the First-Order Approximation (FOA) [1], [20] which consists of a 3D elliptic problem or the Stokes model [6], [15] which consists of a 3D saddle-point problem, the MSSA consists of a system of 2D elliptic equations, and thus is much easier to solve. Moreover, any solver that has been developed for the SSA can be extended to solve the MSSA layer-wise or column-wise. The performance of the resulting numerical solvers is tested for the prognostic benchmark flow-line experiments B and D of the ISMIP-HOM project [22]. In addition, one applies the MSSA model to the first test problem proposed by the Marine Ice Sheet Model Inter-comparison Project (MISMIP) [21].

This paper is organised as follows: in Section 2, the SSA model and its multilayer extension are first recalled. Then the MSSA model is reformulated as a minimisation problem. Afterwards, two numerical methods based on an SSA solver are described for solving the MSSA system in Section 3. Lastly, numerical results are reported in Section 4.