Numerical stability of the free-surface problem arising in ice-sheet modeling

Assessing future sea levels involves computationally expensive simulations of the Earth's glaciers and ice sheets. Computers are, nonetheless, imperfect mathematical machines and inevitably introduce errors into calculations. Central to any numerical algorithm is whether it amplifies or suppresses such errors as the simulation progresses — so-called numerical stability. Stability considerations typically impose restrictions on the time-step size, leading to long simulation times. This is in particular a problem for the most physically accurate and computationally demanding full-Stokes model. In an effort to alleviate stability issues and thereby speed up simulations, this thesis consists of four papers dealing with numerical stability of this system. Papers I and II investigate and adapt a stabilization method from mantle convection to ice-sheet modeling. In these papers, it is shown that the adapted stabilization method substantially increases stable time-step sizes, with an overall negligible impact on accuracy. The results hold for simulations on simpler benchmark domains, as well as for more realistic glacier simulations. The method is easy to implement with minimal impact on computation times, and can therefore be added as a safeguard against stability issues polluting the accuracy of the solver. Paper III is of a more theoretical nature and investigates numerical stability by formulating the coupled system as a single dynamical system, from which numerical stability is analyzed by means of linearization around an equilibrium. Based on this analysis, it is concluded that the increased stability properties of the stabilization follow from a reduction in the spread of the eigenvalues of the linearization. This insight may be used to develop improved stabilization methods based on, for example, state feedback control. In Paper IV, the stabilization is leveraged to develop a higher-order, fully implicit solver — effectively combining stability with accuracy. This is one of the most promising outcomes of this thesis, as stable and accurate numerical methods are necessary for reliable long-term sea-level projections.