  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.
