Effects of numerical implementations of the impenetrability condition on non-linear Stokes flow: applications to ice dynamics

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.

💡 Research Summary

The paper investigates how different numerical implementations of the impenetrability condition at the glacier bed affect solutions of the non‑linear Stokes equations with Navier‑slip boundary conditions, a problem central to ice‑sheet modelling. Three widely used approaches are examined: (1) a “strong” method that locally rotates the coordinate system so that the normal direction aligns with a coordinate axis and then enforces zero normal velocity; (2) a “weak” method that introduces a Lagrange multiplier on each basal facet to constrain the average normal flux; and (3) an “approximative” method that re‑uses the pressure variable—originally the Lagrange multiplier for incompressibility—as a surrogate multiplier for both incompressibility and impenetrability.

The authors first provide a theoretical analysis. The strong method imposes the constraint exactly and does not interfere with the incompressibility condition, but it requires a per‑facet coordinate transformation and associated Jacobian calculations. The weak method adds extra degrees of freedom (the multipliers) to the global system, which can degrade the conditioning of the matrix, especially in three dimensions where the multiplier may become entangled with the pressure field. The approximative method relaxes the incompressibility constraint and enforces impenetrability only insofar as the pressure approximates the normal component of the basal stress (σ·n≈p). Consequently, its accuracy depends on the closeness of pressure to the true normal stress.

Numerical experiments are carried out in two stages. In the first stage, a method of manufactured solutions (MMS) is used to create an artificial exact solution. All three implementations achieve the expected second‑order convergence in L2 and H1 norms, and the differences between them are negligible. This result reflects the fact that the manufactured test does not contain the complex topography and highly variable sliding that characterize realistic ice‑sheet problems.

The second stage employs more realistic two‑ and three‑dimensional test cases with strongly varying basal topography and spatially heterogeneous slip coefficients. In these scenarios the approximative method deviates from the strong method by about 1 % in the 2‑D case and by more than 5 % in the 3‑D case. The error is largest in regions of high basal sliding where the pressure‑stress approximation breaks down. The weak Lagrange‑multiplier method shows acceptable performance in 2‑D but suffers from poor matrix conditioning and slower or failed convergence in 3‑D.

From the combined theoretical and numerical evidence the authors conclude that the strong, locally‑rotated implementation remains the preferred choice for ice‑sheet models. It provides the most reliable enforcement of the impenetrability condition, exhibits stable convergence properties, and integrates cleanly with existing finite‑element Stokes solvers despite the modest overhead of coordinate rotations. The approximative pressure‑based approach, while computationally cheap, should be used with caution: it may be acceptable in regions where pressure is a good proxy for basal normal stress, but it can introduce non‑trivial velocity errors in fast‑sliding or highly undulating beds. The weak method, although conceptually attractive, appears unsuitable for large‑scale three‑dimensional simulations because of its adverse effect on matrix conditioning.

The study thus offers a clear, quantitative guideline for model developers: adopt the strong rotation method as the default, reserve the pressure‑based approximation for exploratory or low‑resolution runs, and avoid the weak multiplier formulation in full‑scale 3‑D applications unless additional stabilization techniques are introduced. Future work could focus on improving the pressure‑based scheme (e.g., by adding corrective terms that better approximate σ·n) or on optimizing the implementation of the strong method to reduce the computational cost of the local rotations.