Constraining the basal sliding of alpine glaciers using adjoint-based inversions and automatic differentiation on GPUs

There are many uncertainties associated with natural processes governing the evolution of glaciers and ice sheets. Unknown parameters, such as basal drag or surface mass balance, can be estimated through data assimilation workflows coupled with physics-based models of ice flow. The multi-scale nature of ice dynamics and significant spatial and temporal variations in physical properties challenges accurate estimations of distributed fields of unknown parameters. High-performance computing and modern computational architectures such as graphics processing units (GPUs) enable efficient and scalable solvers for the ice flow. In this work, we introduce a novel GPU-accelerated inversion framework to enable a point-wise reconstruction of unknown parameter fields at high resolution. The inversion framework is based on the adjoint sensitivity method combined to a gradient-based optimisation. The derivation of the adjoint problem often represents a tedious task which limits the applicability of adjoint-based inversions to simplified ice flow models and hinders fast development. To address these limitations we use differentiable programming and the Julia language which permit automatic differentiation (AD) of arbitrary GPU code. Our GPU-based inversion procedure combines an automatically generated adjoint solver by the Enzyme.jl package using AD and a forward solver to retrieve the point-wise gradient we further use to minimise a cost-function. We demonstrate the capabilities of our inversion framework by developing a forward solver, based on shallow ice approximation (SIA), and several inverse models, utilising different assumptions about the ice flow. One inversion model assumes that the glacier is in a steady-state, which requires iterative solution of SIA equations. The inversion procedure estimates distributions of sliding coefficient, matching the observed ice thickness, glacier outline, and surface velocities. Another inversion model is based on a "snapshot" approach, in which the surface elevation is fixed from observations, and the sliding coefficient is solved to only match observed surface velocities. These two models represent two end members of the spectrum of data assimilation approaches, which will serve as building blocks for more complex workflows, such as transient evolution of glaciers. The feasibility of our inversion algorithms is validated through extensive testing on synthetic glaciers. Then we consider the application of our inversion approach to glaciers in the European Alps using remote sensing data. A map of sliding coefficients is reconstructed by matching the observed ice surface velocity and elevation. The successful application to real glaciers confirms that our inversion models are well suited for large scale and high-resolution simulations. We also present the performance testing results demonstrating close-to-optimal performance of forward and inverse models on NVIDIA GPUs.