  The basal sliding of glaciers and ice sheets can constitute a large part of
the total observed ice velocity, in particular in dynamically active areas. It
is therefore important to accurately represent this process in numerical
models. The condition that the sliding velocity should be tangential to the bed
is realized by imposing an impenetrability condition at the base. We study the,
in glaciological literature used, numerical implementations of the
impenetrability condition for non-linear Stokes flow with Navier's slip on the
boundary. Using the finite element method, we enforce impenetrability by: a
local rotation of the coordinate system (strong method), a Lagrange multiplier
method enforcing zero average flow across each facet (weak method) and an
approximative method that uses the pressure variable as a Lagrange multiplier
for both incompressibility and impenetrability. An analysis of the latter shows
that it relaxes the incompressibility constraint, but enforces impenetrability
approximately if the pressure is close to the normal component of the stress at
the bed. Comparing the methods numerically using a method of manufactured
solutions unexpectedly leads to similar convergence results. However, we find
that, for more realistic cases, in areas of high sliding or varying topography
the velocity field simulated by the approximative method differs from that of
the other methods by $\sim 1\%$ (two-dimensional flow) and $> 5\%$ when
compared to the strong method (three-dimensional flow). In this study the
strong method, which is the most commonly used in numerical ice sheet models,
emerges as the preferred method due to its stable properties (compared to the
weak method in three dimensions) and ability to well enforce the
impenetrability condition.
