Investigation of numerical continuation methods for marine ice-sheet systems formulated as contact problems

Marine ice sheets are complex systems whose response to external forcing is the subject of much attention in the scientific community. In particular, the West Antarctic ice sheet, which could have a significant impact on future sea-level rise, is a major concern. One method of studying the response of marine ice sheets consists in investigating the relationship between the parameters and the equilibrium states of such systems. However, this is typically done by varying these parameters and letting the ice sheets evolve to new steady states, i.e., through transient simulations, which are computationally expensive.

An alternative is to consider continuation methods, where the equilibrium state of a system is studied directly as a function of the parameters. Such an approach has already been used in glaciology to study the mechanical behaviour of 1D marine ice sheets (Mulder et al., 2018), highlighting the hysteresis phenomena that had previously been obtained theoretically (Schoof, 2007). However, this study, because it only considers 1D geometries, does not allow to take into account the effect of lateral drag and of complex bedrock geometries, which are two factors that have the potential to stabilize the grounding line (Gudmundsson, 2013, Sergienko and Wingham, 2021).

Here we consider the continuation problem in the context of 2D marine ice sheets. This introduces several mathematical difficulties, notably related to the treatment of the distinction between grounded and floating parts. Mathematically, the problem takes the form of a contact problem between the bedrock and the lower part of the ice sheet. This leads to a system of equations that is not differentiable, which is challenging to solve numerically. We address these challenges in the context of the continuation problem, and propose several solutions, including a norm-based approach that is inspired from earlier studies (Mittelmann, 1987). Finally, we present some preliminary results which show that our numerical method is promising.