Within the framework of sea-level rise projections, there is a strong need for hindcast validation of the evolution of polar ice sheets in a way that tightly matches observational records (from radar, gravity, and altimetry observations mainly). However, the computational requirements for making hindcast reconstructions possible are severe and rely mainly on the evaluation of the adjoint state of transient ice-flow models. Here, we look at the computation of adjoints in the context of the NASA/JPL/UCI Ice Sheet System Model (ISSM), written in C++ and designed for parallel execution with MPI. We present the adaptations required in the way the software is designed and written, but also generic adaptations in the tools facilitating the adjoint computations. We concentrate on the use of operator overloading coupled with the AdjoinableMPI library to achieve the adjoint computation of the ISSM. We present a comprehensive approach to (1) carry out type changing through the ISSM, hence facilitating operator overloading, (2) bind to external solvers such as MUMPS and GSL-LU, and (3) handle MPI-based parallelism to scale the capability. We demonstrate the success of the approach by computing sensitivities of hindcast metrics such as the misfit to observed records of surface altimetry on the northeastern Greenland Ice Stream, or the misfit to observed records of surface velocities on Upernavik Glacier, central West Greenland. We also provide metrics for the scalability of the approach, and the expected performance. This approach has the potential to enable a new generation of hindcast-validated projections that make full use of the wealth of datasets currently being collected, or already collected, in Greenland and Antarctica.

Constant monitoring of polar ice sheets through remote sensing, in particular
since the advent of altimeter, radar, and gravity sensors such as ICESat-1,
CryoSat, RADARSAT-1, ERS-1 and ERS-2, Envisat, and GRACE, has created a large
amount of data that has yet to find its way through ice sheet models (ISMs)
and hindcast reconstructions of polar ice sheet evolution. In particular, as
evidenced by the wide discrepancy between ISMs involved in the SeaRISE and
Ice2Sea projects

Applications to transient models and long temporal time series such as the
ICESat/CryoSat continuous altimetry record from 2003 to the present day have
been much more rare, and, to our knowledge, are limited to a few studies such
as

Understanding sensitivities of a forward ice-flow model, which is needed to
physically constrain a temporal inversion, requires computation of
derivatives of model outputs to model inputs. If such derivatives are
approximated by finite-difference schemes, they are subject to the tradeoff
between approximation and truncation errors for the perturbation, which is
aggravated for higher-order derivatives. If the derivatives are computed
using algorithmic differentiation (AD;

In practice most

This is achieved using the partial derivatives of the forward operation

For the multiplication example one therefore gets

Using this approach, it is possible, without manual derivation, to evaluate
the derivative

Applying AD to large-scale frameworks such as the Ice Sheet System Model
(ISSM)

Distinctions between the AD implementation approaches relate to the way the
derivative data in

Changing to the aforementioned active type is a significant effort to be undertaken in the model code. Among the choices to effect this type change, one should select one that is transparent to the model development process, is maintainable, and minimizes the manual effort.

The tool of choice for this paper is ADOL-C

Instead, the recommended approach is the use of C++ template classes and
methods where

Overall, the use of C++ template classes introduces an additional element of
complexity to the code, but it also means that the same exact framework can
be used to compute both the forward run and its gradient, using the exact
same code base, without the need to manage different code for the computation
of

Interpreting the trace by

Thus, the use of

The strength of the AD approach lies in the automated differentiation of all
the parts of the model that are changing and expanding as the model is being
developed, as opposed to manually programming adjoints that are error-prone,
and implies a significant effort in terms of code development and
maintenance. However, because models reference libraries, for fixed
mathematical mappings that may have easily programmable high-level adjoints,
one can and should exploit this insight as it can yield efficiency gains not
obtainable by an AD tool remaining unaware that calls to some collection of
subroutines represent a certain mathematical mapping. A good example are
linear solver libraries that may not even be written in the same programming
language as the model in question. Considering the system

Is

Are the previous values of

If

How does the external function interface used by

External functions

expanded the

added optional parameters to pass in the sparsity pattern of

added controls for storing and restoring prior values of

added tracking of the maximal

Regarding Q1, we know that in the ISSM,

Logic variants encapsulated in the ISSM MPI wrapper library.

As is the case with the ISSM, practically relevant science problems incur a
computational complexity that necessitates execution on parallel hardware,
often using MPI as the paradigm of choice. Sending data with MPI from a
source buffer

For testing and small-scale experiments, the ISSM code, as many other models,
treats MPI parallelization as a compile-time option controlled by
preprocessor macros. Furthermore, turning the adjoint capability on and off
as suggested in Sect.

The approach is shown for

Code for wrapping reduction.

Code snippets for templating with corresponding MPI data types.

Variant 1 is implemented by line 21, variant 2 by lines 13–18 and 21,
variant 3 by line 9, and variant 4 by line 7. The example code reflects some
of the simplifying assumptions being made based on the MPI usage patterns in
the ISSM. Examples are the absence of

Types per variant from Table

ISSM AD validation suite integrated within Jenkins for continuous integration and delivery

The matching between the actual type of the buffer and the corresponding MPI
type must extend to the templating approach suggested for the type change in
Sect.

Gradient of the surface altimetry cost function (spatio-temporal
average of the misfit between the 6-year record of observed ICESat altimetry
and modeled surface elevation) with respect to SMB (in years). The gradient
is computed for every time step of a 2-week interval time series between 2003
and 2009. We only show three periods (September 2003, June 2006, and February
2009) in the entire interval. The location for the study is the northeastern
Greenland Ice Stream, and the gradient is taken from an inversion study of
the surface forcings necessary to best fit the ICESat altimetry record

Gradient of the surface velocity cost function (spatio-temporal
average of the misfit between a 5-year time series of observed surface
velocities and the modeled surface velocity) with respect to the basal drag
coefficient (in yr

Finally, a practical concern for using the parallelized adjoint is the
handling of sensitivities to quantities that are uniformly initialized across
ranks, such as parameters. Frequently, as was the case within the ISSM, these
quantities are initialized from files or otherwise per process in the
parallel case the same way as in the sequential case. In the parallel case
that implies a replication of the same quantity across ranks. However, to
obtain the correct sensitivities, the quantity

The ISSM is validated in AD mode by continuously running a test suite within
the Jenkins

Test problem setup with the sequential GSL solver; plots over max mesh distance show overhead factors (upper frame) and approximate runtime portions of the model execution with adjoint computation (lower frame).

Test problem setup with the parallel MUMPS solver; plots over max mesh distance show approximate runtime portions of the model execution with adjoint computation (upper frame) and overhead factors for the trace creation and interpretation (lower frame).

The ongoing use of adjoint computations includes sensitivity studies and
state estimation problems for transient model runs. Because this paper
concentrates on technical aspects, we show, only as exemplary evidence of the
practical usefulness, some sensitivities of cost functions calibrated for two
sensors commonly encountered in cryosphere science, altimeters (that measure
surface elevation) and radars (that measure surface displacement or
velocity). In Fig.

A second type of observation garnering considerable interest is temporally
resolved time series of radar observations (using speckle-tracking or InSAR
to infer surface deformation of the ice) to measure variations of surface
velocity on a seasonal timescale (one observation every 2 months, as in

Currently, these adjoint computations are performed on the Pleiades cluster
at the NASA Advanced Supercomputing Center. Given the fact that the effort
for computing the adjoint

The earliest ISSM adjoint computations took place before the MPI wrapper
library was introduced and therefore were done sequentially with the GSL
solvers. To evaluate the performance, a representative test case was chosen
from the ISSM regression test suite, test 101. This test models the
steady-state stress balance of a square ice shelf of size
50 km

In combination with the much reduced impact of the solver, the overhead
factor for equivalent test problems reached temporarily up to 145, which
clearly was not acceptable. Subsequent changes to the ADOL-C tool included
improvements in the internal address management. Changes to the model
included modifications of the sparse data format, enabling the control of the
I/O buffering for the trace

We developed a new adjoint capability in the ISSM, based on the ADOL-C framework and the AdjoinableMPI library, which is to our knowledge the first time this type of approach has been systematically applied to a software framework of this size and complexity. Despite the difficulties encountered in rewriting the software, the overloaded approach is transparent to the user, which is critical given the size of the larger cryosphere science community that is not familiar with the adjoint work, and for which classic approaches such as source-to-source transformation may not be applicable or may prove more cumbersome. The flexibility of this approach allows in particular for quick turn-around in developing adjoint models of new parameterizations that are not easily hand-derived. This is a major advantage in that it opens this approach to the wider community. This, given the large amount of remote sensing data currently being collected and under-utilized, could prove paramount if we are to hindcast validate projections of sea-level rise. Further work is of course required to bring in additional observations such as gravity sensors, or radar stratigraphy observations, which will involve development of new cost functions, and scalability in 3-D. Though this is complex in that it requires integrated resiliency and adjoint checkpointing schemes for long-running transient modeling scenarios, our approach has proven flexible, and should lead to a brand new set of data assimilation capabilities that have already been available to other Earth science communities for a long time. Indeed, by allowing temporal data assimilation for a large number of sensors and models, such as demonstrated here with the use of altimetry and radar sensors for mass transport and stress balance models, respectively, the ISSM paves the way for wider integration between the modeling and observational cryosphere communities.

The ISSM code and its AD components are available at

DE-AC02-06CH11357. Eric Larour and Jean Utke were supported by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the NASA Cryospheric Sciences and Modeling and Analysis Programs. Gilberto Perez was supported by a subcontract from the Jet Propulsion Laboratory to the University of California at Irvine and Mathieu Morlighem was supported under a contract with the NASA IceBridge research program. Edited by: D. Goldberg Reviewed by: L. Hascoet and one anonymous referee