Comparing methods of computing second order derivatives in numerical ice sheet models

Large-scale continuum ice sheet models solve a set of conservation-of-momentum equations to calculate ice speed given the ice sheet geometry and, often, some a priori unknown parameters. To facilitate the initialisation of projections, many such models allow one to calcuate gradients of model outputs with respect to those unknown input paramters. In some models, this is done with the help of algorithmic differentiation (AD), while in others it is done using hand-derived PDE-level adjoint rules. Increasingly, there is interest in computing higher-order derivatives, for example, to facilitate the use of different optimisation algorithms or perform uncertainty quantification. In this work, we derive and implement a second-order adjoint (SOA) model for the shallow stream equations, implementable in any 2D ice sheet model. We implement this in a finite volume code written in a numerical computation library for Python called JAX. We conduct comparisons between the computation of Hessians using this SOA model and using the AD tools provided by JAX. The SOA model makes the assumption of linear rheology, as many first-order adjoint models do. We find that this causes a rapid departure in the directions of the principal components of the Hessian from those founding using AD. Hence, we consider AD to be the more suitable choice for many applications.