Galerkin Least-Squares Stabilization in Ice Sheet Modeling - Accuracy, Robustness, and Comparison to other Techniques

TL;DR

This paper evaluates the accuracy and robustness of Galerkin Least-Squares stabilization in ice sheet modeling, compares it with other methods, and discusses its sensitivity and limitations near ice margins.

Contribution

It provides a comprehensive comparison of GLS stabilization with other techniques and highlights its limitations and sensitivities in ice sheet simulations.

Findings

Vertical velocity is more sensitive to stabilization parameters.

Standard GLS parameters may cause oscillations near ice margins.

Interior penalty method offers better accuracy in challenging regions.

Abstract

We investigate the accuracy and robustness of one of the most common methods used in glaciology for the discretization of the $\mathfrak{p}$-Stokes equations: equal order finite elements with Galerkin Least-Squares (GLS) stabilization. Furthermore we compare the results to other stabilized methods. We find that the vertical velocity component is more sensitive to the choice of GLS stabilization parameter than horizontal velocity. Additionally, the accuracy of the vertical velocity component is especially important since errors in this component can cause ice surface instabilities and propagate into future ice volume predictions. If the element cell size is set to the minimum edge length and the stabilization parameter is allowed to vary non-linearly with viscosity, the GLS stabilization parameter found in literature is a good choice on simple domains. However, near ice margins the standard parameter choice may result in significant oscillations in the vertical component of the surface velocity. For these cases, other stabilization techniques, such as the interior penalty method, result in better accuracy and are less sensitive to the choice of the stabilization parameter. During this work we also discovered that the manufactured solutions often used to evaluate errors in glaciology are not reliable due to high artificial surface forces at singularities. We perform our numerical experiments in both FEniCS and Elmer/Ice.