SICOPOLIS-AD v2: tangent linear and adjoint modeling framework for ice sheet modeling enabled by automatic differentiation tool Tapenade

SImulation COde for POLythermal Ice Sheets ( SICOPOLIS ) is an open-source, 3D dynamic/thermodynamic model that simulates the evolution of large ice sheets and ice caps. SICOPOLIS has been developed continuously and applied to problems of past, present, and future glaciation of Greenland, Antarctica


Statement of need
The two contemporary ice sheets, Greenland and Antarctica, are dynamic entities whose evolution is governed by a set of nonlinear partial differential equations (PDEs) that describe the conservation of mass, momentum, and energy, as well as constitutive laws for the material properties of ice. In general, these equations cannot be solved analytically but must be solved numerically. Ice sheet models are a computer representation of these PDEs. They require as input parameters (i) initial conditions of the state of the ice sheet, (ii) surface boundary conditions, such as precipitation, (iii) basal boundary conditions, such as geothermal flux, and (iv) model parameters, such as flow law parameters. Despite advances in numerical modeling of ice sheets, the effects of ad-hoc initialization and the uncertainties in these independent input parameters propagate to quantities of interest (QoI), such as future projections of sea-level rise, which is of economic and societal importance (Schinko et al., 2020). It is thus desirable to evaluate the sensitivities of our QoI to these independent input variables.
In the context of ice sheet modeling, sensitivities of model-data misfits or other QoI are a key ingredient for performing model calibration, state estimation, or uncertainty quantification (UQ), which guide the improvement of model simulations through PDE-constrained gradient-based optimization.  leverages the recently open-sourced AD tool Tapenade (Laurent Hascoët & Pascual, 2013) to generate code for the adjoint model of the open-source ice sheet model, SICOPOLIS (Greve, 1997;Greve et al., 2011;Greve & Blatter, 2009). Sensitivities can be calculated using a single forward and adjoint model evaluation, instead of the ( ) forward model evaluations. Empirically, one adjoint model evaluation is about 5-10 times as expensive as a forward model run. The adjoint computation is highly efficient for calculating sensitivities when is large (typically, ∼ 10 4 − 10 6 ).
The functionality to generate a tangent linear version of the forward model is also included, which was not available in SICOPOLIS-AD v1. This is valuable for UQ of the inferred parameters, as well as uncertainty propagation to QoIs. It can also be used to verify the results of the adjoint model.

Target Audience
This package is intended as a resource that enables sensitivity analysis, model calibration, and uncertainty quantification of a continent-scale ice sheet model. Our package is also intended to serve as a guide for future work in the application of open-source AD tools for physics-based simulation codes written in Fortran.

State of the field
SICOPOLIS is among the early thermo-mechanical models to simulate contemporary and paleo continental-scale ice sheets (Greve, 1997). Like similar models developed at the time, including Glimmer and its successor, the Community Ice Sheet Model (CISM) (Rutt et al., 2009), GRISLI (Ritz et al., 1996), the model by Huybrechts (1990), or by Pollard & DeConto (2009), SICOPOLIS has been based (until recently) on the so-called shallow ice approximation to simplify the Cauchy stress tensor in the momentum conservation equation, implemented on a regular, finite-difference mesh. See Hindmarsh (2004) for other approximations commonly employed in ice sheet models. This approximation enabled the efficient computation of ice sheet evolution over long, glacial/deglacial cycles.
The last decade has seen substantial advances in continental-scale ice sheet modeling, with the development of several new ice sheet models (some of which are on unstructured grids using finite element methods), notably the Ice Sheet System Model ISSM (Larour et al., 2012), the Parallel Ice Sheet Model PISM (Bueler et al., 2007), Elmer/Ice (Gagliardini et al., 2013), or the MPAS-Albany Land Ice MALI (Hoffman et al., 2018). While designed to capture the evolution of short-term, fast-flowing, or fast-changing outlet glaciers via horizontal stress contributions, these models have so far found little application in paleo-ice sheet simulations due to their extensive computational costs. A compilation of the suite of ice sheet models used for the latest Ice Sheet Model Intercomparison Project, Phase 6 (ISMIP6) in support of the IPCC's Sixth Assessment Report is available in Payne et al. (2021) and Nowicki et al. (2016).
Relevant to this paper, of all the time-evolving models listed, apart from SICOPOLIS-AD (Heimbach & Bugnion, 2009;Logan et al., 2020), only the ISSM model and variants thereof possess adjoint model codes which have been generated, in part, using automatic differentiation (L. Hascoët & Morlighem, 2018;Larour et al., 2014). Multi-centennial and longer integrations with the adjoint model have so far been conducted only with SICOPOLIS-AD.

Features
AD tools such as the commercial TAF (Giering & Kaminski, 1999) and the open-source OpenAD (Utke et al., 2008) have been used previously with SICOPOLIS (Heimbach & Bugnion, 2009;Logan et al., 2020). OpenAD is no longer actively developed because it is based on the Open64 compiler which ceased development in 2011. The differentiation of SICOPOLIS, therefore, must be performed using a different tool. Compared to OpenAD, the Tapenade enabled implementation has the following advantages: • It is up-to-date with the latest SICOPOLIS code.
• The AD tool Tapenade is open-source and actively maintained.
• A new tangent linear code generation capability is introduced.
• The AD-generated codes can accept NetCDF inputs.
• The external library LIS, its tangent linear code, and adjoint code are correctly incorporated which can improve the simulation of Antarctic ice shelves and Greenland outlet glaciers.
• Gitlab-CI, a Docker, and the pytest framework are leveraged for Continuous Integration (CI) to track changes in the trunk that "break" the AD-based code generation.
• The entire code is parsed by Tapenade, preventing cumbersome manual maintenance of subroutines to initialize the adjoint runs.
• Python scripts are provided for quick setup of the compilation, I/O, and execution processes based on user-provided metadata.
• The setup is well-documented, along with tutorials.

Software requirements and external usage
SICOPOLIS-AD v2 is built on top of the ice sheet model SICOPOLIS and uses Tapenade to differentiate this model. All the prerequisites of using SICOPOLIS and Tapenade need to be satisfied. A Python installation is needed to use the automation tools.

Example
We illustrate the use of our tool with the example of a steady-state simulation of the Greenland ice sheet under modern climate conditions. The corresponding SICOPOLIS configuration header file, v5_grl16_bm5_ss25ka, is provided as a reference template in the standard SICOPOLIS distribution. We shorten the total integration time to 100 simulated years to keep the computational cost of the tangent linear and finite differences reasonable. Our QoI (i.e., dependent variable) is the total volume of the ice sheet at the end of the run (fc). The sensitivity is evaluated with respect to the geothermal heat flux, q_geo (independent variable), a 19,186-dimensional field. The results are shown in Figure 1. The results show good agreement between all three modes used to evaluate this sensitivity. The error is less than 6% between AD-generated (adjoint/tangent linear codes) and finite differences at all but one point with "significant" gradient values, i.e. within 4 orders of magnitude of the maximum absolute value of the finite differences gradient. The relative error between the AD-generated adjoint and tangent linear models is less than 0.002% at all points with values within 4 orders of magnitude of the maximum absolute value of the finite differences gradient. However, the adjoint model is much faster than the other two, as shown in Table 1, because the number of evaluations of the latter two scales linearly with the parameter dimension (~( )). The discrepancy will be even larger if a finer mesh is used.