Generally, we consider maps from a set of parameters to a set of quantities of interest that we seek to calibrate to observational data. In particular, we focus on maps that incorporate physical models in form of PDEs. Computing the solution of the so-called parameter-to-observable map constitutes a forward problem, while finding parameters such that the model output is consistent with observational data amounts to an inverse problem.

As an application, we target Earth's mantle convection with associated plate tectonics at a global scale. The flow of (solid) rock in the mantle is modeled by nonlinear, incompressible Stokes equations. These complex models can only be solved numerically and are computationally expensive requiring: (a) advanced numerical methods with (nearly) optimal algorithmic scalability (with increasing degrees of freedom) and (b) parallel scalability (with increasing number of compute cores). Given a set of model parameters, the (forward) simulation of Earth's mantle convection produces a global velocity field. These velocities are compared to observed plate velocities at Earth's surface. The inverse problem consists of finding a set of parameters whose model output at the surface is consistent with the observational data.

We pose the inverse problem in a Bayesian framework by introducing a prior for the parameters. We compute the maximum a posteriori (MAP) estimate and an approximation of the parameter uncertainties. This amounts to the solution of an optimization problem governed by the model equations, and computation of second derivatives at this MAP point. Newton’s method is used for solving the optimization problem, which requires first- and second-order derivatives (gradients and Hessians, respectively). These operations are performed in an efficient and scalable way using adjoint methods.