## Inverse problems for basal properties in a thermomechanically coupled ice sheet model

##### Abstract

Polar ice sheets are losing mass at a growing rate, and are likely to make a dominant contribution to 21st century sea-level rise. Thus, modeling the dynamics of polar ice sheets is critical for projections of future sea level rise. Yet, this is difficult due to the complexity of accurately modeling ice sheet dynamics and because of the unobservable boundary conditions at the base of the ice sheet.
In this work, we address the inverse problem of inferring basal properties--in particular, the geothermal heat flux--from surface velocity observations and a forward model in the form of thermomechanically coupled nonlinear Stokes and energy equations with complementarity conditions representing melting of basal ice. This inverse problem is made even more challenging due to the severely nonlinear and non-smooth nature of the forward problem. The inverse problem is formulated as a nonlinear least squares optimization that minimizes the misfit between the model prediction and the observation. A Tikhonov regularization term is added to render the problem well-posed. To solve the inverse problem for large-scale ice sheets, the use of adjoint-based Newton methods is critical, which requires a smoothly differentiable forward problem. Thus, we regularize the complementarity conditions of the forward problem by a penalty-like method, such that the solution of the regularized problem approaches that of the original forward problem as the penalty approaches infinity. As a consequence of the Petrov-Galerkin discretization of the energy equation, discretization and differentiation do not commute, i.e., the order in which we discretize the cost functional and differentiate it affects the correctness of the gradient. Here, we employ the discretize-then-optimize approach to guarantee the consistency between the discrete cost function and its gradient. Using two- and three-dimensional model problems, we study the prospects for and limitations of the inference of the geothermal heat flux field from surface velocity observations. The results show that we can approximately locate the melting region of the ice sheet and reconstruct the geothermal heat flux in the cold region. The reconstruction improves as the noise level in the observations decreases but short-wavelength variations in the geothermal heat flux are difficult to recover. We analyze the ill-posedness of the inverse problem as a function of the number of observations by examining the spectrum of the Hessian of the cost functional. Motivated by the popularity of operator-split or staggered solvers for forward multiphysics problems---i.e., those that drop two-way coupling terms to yield a one-way coupled forward Jacobian---we study the effect on the inversion of a one-way coupling of the adjoint energy and Stokes equations. We show that taking such a one-way coupled approach for the adjoint equations leads to an incorrect gradient and premature termination of optimization iterations. This is due to loss of a descent direction stemming from inconsistency of the gradient with the contours of the cost
functional.
For large-scale simulations, we extend the parallel solver "ymir" for the nonlinear Stokes model by incorporating the coupled nonlinear advection-diffusion energy equation with complementarity conditions. An inexact Newton-Krylov method with block preconditioning is designed for the thermomechanically coupled ice sheet model. The inexact Newton-Krylov method exhibits near optimal algorithmic scalability, i.e., the numbers of both Newton and Krylov iterations depend only weakly on problem size. We use the parallel solver to simulate the flow of the Antarctic ice sheet and the Pine Island ice stream, with a geometry and boundary conditions from the ALBMAP dataset.