Numerical solution of a thermomechanical coupled model governing glacier evolution

In this paper the numerical solution of a highly nonlinear model for the thermomechanical behavior of polythermal glaciers is presented. The modeling follows the shallow ice approximation (SIA) for glaciers introduced in Fowler (1997)  [13]. The model has been extended to incorporate additional moving boundaries and other nonlinear features. Moreover, a fixed domain formulation is proposed to avoid the computational drawbacks of a time-dependent domain in the numerical simulation with front tracking methods. In this setting, the coupled problem is decomposed into different nonlinear problems which allow one to obtain sequentially the profile evolution, the velocity field, the glacier surface and atmospheric temperatures, basal magnitudes and the temperature distribution inside the ice mass. A fixed point iteration algorithm converges to the solution of the nonlinear coupled problem. Among different numerical methods involved in the solution of the subproblems, characteristic schemes for time discretization, finite elements for spatial discretization, duality methods for the nonlinearities associated to maximal monotone operators and a Newton scheme for the nonlinear viscous term are proposed. Several numerical simulation examples illustrate the performance of the numerical methods and the behavior of the involved physical magnitudes.