Numerical ice sheet modelling using FEniCS

The increase in computational power in the last 20 years has made it feasible to solve the full Stokes equation for ice flow, using finite element methods (FEM). However, the numerical properties of these equations remains largely unknown. This is due to their nonlinear nature, which makes them hard to analyze mathematically. For this reason convergence rate and stability has to be established by performing simulations. In this thesis a 2D ice solver has been developed. The solver was then tested on the ISMIP-HOM benchmarks, in order to assert convergence rates and stability. The solver was developed using the FEM software FEniCS. It was verified that the solver could obtain the same convergence rate as in the case of linear Stokes; however, for curved boundaries the convergence rate in velocity dropped by one order. It was also found that discontinuity in the coefficient of basal traction reduced convergence rate to linear in both pressure and velocity. Time dependence was added to the solver by coupling the Stoke's system to the kinematic free surface equation. Time dependence was then tested using one of the EISMINT benchmark. It was found that for a mesh with a mesh size parameter of 10 km, instabilities arised after 2500 years when using a time step of 25 years, resulting in spurious numerical oscillation. However, by introducing a diffusive term into the surface equation, it was possible to attenuate these oscillation without affecting the overall shape of the ice sheet