2.1 Forward model SICOPOLIS

We begin with the ice sheet model SICOPOLIS (SImulation COde for POLythermal Ice Sheets) and sketch the development of its adjoint model from version 5-dev (Greve, 2019; Rückamp et al., 2019) (http://www.sicopolis.net, last access: 2 April 2020). SICOPOLIS is open source and written in Fortran; it has a relatively long and stable history (Greve, 1997). SICOPOLIS has remained a relevant and powerful tool for the cryosphere community and continues to participate in model intercomparison exercises (e.g., Goelzer et al., 2018; Seroussi et al., 2019). The model couples ice sheet dynamics and thermodynamics (solving for the ice thickness, extent, velocity, temperature, water content, and age) over three-dimensional domains including, among others, Greenland, Antarctica, paleo-ice sheets on Earth, and the polar caps of Mars. It employs simplified versions of the three-dimensional Stokes equation (internal stress balance). These include the shallow-ice approximation for ice resting on land (Hutter, 1983; Morland, 1984), the shallow-shelf approximation for ice floating in the ocean (Morland, 1987; MacAyeal, 1989; Weis et al., 1999), and the shelfy-stream approximation for fast-flowing ice streams with limited coupling to the bed (Bernales et al., 2017). For a detailed treatment of the numerical methods employed in the model, readers are referred to Greve and Calov (2002), Greve and Blatter (2009, 2016), and Bernales et al. (2017).

SICOPOLIS employs four different thermodynamics representations: (1) a two-layer polythermal scheme, which allows for the computation and effects of liquid water within a warmer temperate layer; (2) a purely cold-ice scheme in which no liquid water is present; and (3, 4) two flavors of the one-layer enthalpy scheme that combine the physical adequateness of (1) with the greater numerical simplicity of (2) (Aschwanden et al., 2012; Greve and Blatter, 2016). In all cases, horizontal diffusion of the thermodynamic fields (temperature, water content, or enthalpy) is neglected. The solvers employed use an implicit discretization scheme for the vertical derivatives and an explicit scheme for the horizontal derivatives.

SICOPOLIS simulates ice as a nonlinear viscous fluid by employing Glen's flow law (Glen, 1955) amended as in Greve and Blatter (2009):

where η is the ice viscosity, T^′ is the temperature difference relative to the pressure-melting point, σ_e is the effective shear stress, and σ₀ is a small constant used to prevent singularities when σ_e is very small. n is the flow law exponent (taken as 3), and A is a temperature- and pressure-dependent rate factor (Cuffey and Paterson, 2010) that is modified in temperate regions containing liquid water following Lliboutry and Duval (1985).

Basal sliding under grounded ice links the sliding velocity, v_b, to the basal shear traction, τ_b, and the basal normal stress, N_b (counted positive for compression), in the form of a Weertman–Budd-type sliding law (e.g., Weertman, 1964; Budd et al., 1984):

where C_b is the sliding coefficient, and p and q are the sliding law exponents.

2.2 Adjoint model of SICOPOLIS generated with OpenAD

As described in Sect. 1.1, the construction of an adjoint model of a nonlinear, time-dependent forward model often presents a formidable task when solved analytically or hand-coded (e.g., Goldberg and Sergienko, 2011; Gillet-Chaulet et al., 2012; Morlighem et al., 2013; Isaac et al., 2015). In previous works, the variational forward and adjoint equations are derived first and then discretized. As an alternative, AD produces adjoint code through the differentiation of source code using source transformation or operator overloading tools. As the standard of numerical models (in various contexts) has risen to more complicated physical representations, the use of AD has become increasingly popular (e.g., Heimbach and Bugnion, 2009; Goldberg et al., 2016; Hascoët and Morlighem, 2018, using source transformation – Brinkerhoff and Johnson, 2013; Larour et al., 2014; Hoffman et al., 2018, using operator-overloading and code-templating).

The adjoint of the ice sheet model SICOPOLIS is largely generated automatically by the application of the freely available source transformation tool OpenAD (Utke et al., 2008) developed at Argonne National Laboratory, the University of Chicago, and Rice University (http://www.mcs.anl.gov/OpenAD, last access: 2 April 2020). It is a flexible and modular tool that parses a given model written in Fortran to generate a Fortran version of the model's adjoint code. As opposed to adjoint models that are analytically derived and discretized, the adjoint model of SICOPOLIS furnished by AD yields not one unique adjoint code but many possible configurations of the adjoint model, whose constituent parts are selected at compile time by options in the header files. Thus, there no single, core adjoint code; there is instead a flexible adjoint code that evolves in tandem with the forward model of SICOPOLIS and can accommodate updated empirical relationships and boundary conditions.

The adjoint model of SICOPOLIS produced by OpenAD results in approximately 50 000 executable lines, represented in a much simplified schematic in Fig. 1 by the composition of the blue, red, and black algorithmic steps. Like many complex, time-evolving geophysical models, SICOPOLIS comes with a range of choices of model configuration, in particular numerical schemes, which the user may choose from. As a matter of convenience, the preferred implementation is to make all of these choices (or options) available at runtime to minimize the need for recompiling the model. The same convenience is available, in principle, to the AD-generated adjoint model. The control flow analysis of the AD tool identifies all possible flows of forward model execution and produces corresponding adjoint flow paths. However, close to 2 decades of experience with the application of AD to complex, time-evolving geophysical models, all of which have a range of numerical schemes that users may choose from (Heimbach et al., 2002, 2005; Forget et al., 2015), has shown that for the specific application of adjoint modeling, it is preferable to remove code that will not be executed in a given application from adjoint code generation (and subsequent compilation). The two main reasons for proceeding in this manner are the following.

• i.

  The exclusion of forward model code that the user knows will not be executed may significantly simplify the AD tool's dependency and flow control analysis, avoid spurious dependencies that the AD tool may detect, and lead to more streamlined source code for the adjoint.

• ii.

  Because of the reverse mode and requirement to store required variables in time-reversed order (e.g., those used for evaluating state-dependent conditions and nonlinear expressions), adjoint models will have a substantially larger memory footprint than their parent forward model (Heimbach et al., 2005). Memory requirements may be significantly increased if the adjoint model is required to keep track of a large range of conditional branches for execution.

For these practical considerations, removing nonused forward model code at the time of adjoint code generation and subsequent compilation has proven to be highly preferable (although not strictly required). It is implemented here via C preprocessor (CPP) options that are enabled or disabled prior to generating the adjoint code (and prior to compilation time). We note that the implementation keeps runtime parameters and flags in place such that the forward model default to keep all code available at runtime is not compromised. By pairing SICOPOLIS with the source transformation tool OpenAD, the adjoint model of SICOPOLIS may be generated automatically for a large variety of forward model configurations (including detailed choices of model domain, numerics, control variables, and QoI).

A number of algorithmic aspects of the code needed one-time editing or refactoring for OpenAD to be able to successfully parse the source code and provide correct adjoint code. For example, non-smooth functions – such as piecewise linear functions represented by if statements or absolute values – are inherently non-differentiable and sometimes required special treatment before the adjoint could be obtained by AD. Because of its importance in the development of a forward model that works properly within the framework of the AD tool, we have devoted a detailed description in Appendix B of the aspects of SICOPOLIS that required code refactoring. Further technical details on how to set up, compile, and run reference configurations are documented in a quick-start manual (Logan et al., 2019).

3.1 Antarctic model configuration

We simulate Antarctica for 100 years of model time with a 20 km horizontal resolution and 81 terrain-following vertical layers. The dynamic and thermodynamic time steps (which can be chosen to differ) were both set to 0.2 years, as this was found to be the most stable value for the forward model simulation. Land ice, floating ice, and ice streams are approximated by the shallow-ice approximation (SIA), shallow-shelf approximation (SSA), and shelfy-stream approximation (SStA) formulations, respectively, described in Bernales et al. (2017). Ice thickness evolves freely and without adjustment. Solutions to the SSA portion of SICOPOLIS are aided by invoking the external Library of Iterative Solvers (LIS; https://www.ssisc.org/lis/, last access: 2 April 2020). Thermodynamics are formulated by the conventional enthalpy scheme (Sect. 2.1). We use Glen's flow law (Eq. 8) with a stress exponent n=3, a residual stress σ₀=10⁴ Pa, uniform flow enhancement factors E=5 for grounded ice and E=1 for floating ice, and a rate factor A(T^′) as in Cuffey and Paterson (2010). Horizontal and vertical advection in the temperature and age equations are discretized via a first-order upstream stencil of interpolated velocities and advection terms on the main grid, and topography gradients are evaluated with a fourth-order discretization. The ice temperature is initialized as a uniform value of −10 ^∘C, as the goal of this exercise is a proof of concept and not an exhaustive examination of all aspects of the Antarctic Ice Sheet. For the same reason, a uniform geothermal flux of 55 m W m⁻² is applied. Parameterization of the mean annual and mean January surface temperatures is according to Fortuin and Oerlemans (1990), and the applied surface temperature is held constant throughout the simulation. Accumulation is applied at present-day levels throughout the simulation (Le Brocq et al., 2010; Arthern et al., 2006). The fraction of solid precipitation is a linear function of the monthly mean surface temperature according to Marsiat (1994). Surface ablation is parameterized by the positive degree day method, and rainfall is assumed to run off instantly. Floating ice is removed at calving fronts for thicknesses less than 30 m. The parameters for the basal sliding law (Eq. 9) are chosen as C_b=11.2 m yr⁻¹ Pa⁻¹, p=3, and q=2. Basal melting under floating ice is set to 30 m water equivalent per year around the grounding zone (adjacent grounded and floating grid points) and zero elsewhere, for simplicity. Sea level is constant, and there is no special treatment of subglacial hydrology. The initial geometry is taken from the present day (Fretwell et al., 2013). No isostatic bedrock adjustment is simulated in our configuration. No special pretreatment or “spin-up” of the model is applied prior to the 100-year simulation; rather, the above conditions are set and the model evolves for 100 years.

We note that the operation underlying calving (see above) amounts to a conditional statement. From an AD perspective, the following steps occur: (i) derivative codes are generated for each condition; (ii) code to store and restore the required variable is added to properly evaluate the conditional derivative. For legacy code the operation may not be differentiable at the exact condition (see Appendix B for practical details). This should be taken into account when performing gradient-based optimization.

3.2 Results

The motivation for developing an adjoint of a numerical model stretches far beyond providing comprehensive sensitivity experiments; often, an adjoint model is developed so that the sensitivities may be used in a gradient-based PDE-constrained optimization problem to invert for uncertain initial conditions, boundary conditions, or model parameters, thereby producing a data-constrained estimate for the evolution of the state of the system. Here, however, we are interesting in understanding model sensitivities. We present the sensitivity of the volume of the Antarctic Ice Sheet with respect to several control variables as a proof of concept, rather than extending the work in the direction of optimization, which will be the subject of future studies. The purpose is to gain physical insight into the model's linear response characteristics and to ascertain correctness and interpretability of the adjoint. The adjoint-derived sensitivities are compared to finite-difference perturbations, either at single points or over a patch of the domain that has been uniformly perturbed, to demonstrate that the adjoint model is sufficiently consistent with sensitivities derived via finite differences. Those comparisons are shown in Table 1. Lastly, compared to Heimbach and Bugnion (2009), we present this work for the novel application of examining Antarctic-wide adjoint-generated sensitivity maps. To the authors' knowledge, such a presentation has not been formally examined heretofore.

Figure 2Adjoint sensitivities, $\frac{\mathit{\delta }J}{\mathit{\delta }{X}_i}$, for the Antarctic Ice Sheet, where X_i is the control variable shown. Control variables are the (a) initial ice thickness (m²), (b) mean January precipitation (m² yr), (c) surface temperature (m^3 ∘C⁻¹), and (d) basal temperature (m^3 ∘C⁻¹). Locations in (a) numbered 1–4 are compared to finite-difference values in Table 1.

Figure 3Logarithms of the absolute value of adjoint sensitivities, $\log \left|\frac{\mathit{\delta }J}{\mathit{\delta }{X}_i}\right|$, for the Antarctic Ice Sheet, where X_i is the control variable shown. Control variables are the (a) initial ice thickness (m²), (b) mean January precipitation (m² yr), (c) surface temperature (m^3 ∘C⁻¹), and (d) basal temperature (m^3 ∘C⁻¹).

Figures 2 and 3 respectively show the raw and logarithm of the absolute value (${\log }_{\mathrm{10}}|•|$) of the adjoint sensitivities. We have chosen to present the adjoint values in both ways so that the general pattern and sign of the adjoint values are readily apparent (Fig. 2) as is the order of magnitude of the adjoint values (Fig. 3), which can vary widely across the Antarctic Ice Sheet depending on the control variable. Further, we have chosen the locations shown in Fig. 2a so that the included dynamics and solvers invoked in the code can be tested in three different and important regions and regimes: location 1, the fast-moving Thwaites Glacier, which directly discharges into the Amundsen Sea Embayment; location 2, the middle of the Ross Ice Shelf; and location 3, Slessor Ice Stream, which feeds the Ronne–Filchner Ice Shelf. Testing the agreement between adjoint and finite-difference values at these locations offers a broad sense of the performance of the adjoint model in very different and important environments and dynamical regimes across the Antarctic Ice Sheet. The control variables we have selected to test involve either initial (ice thickness) or boundary conditions (summer precipitation, surface and basal temperature) and are independent inputs to either the conservation of mass (ice thickness and summer precipitation) or conservation of enthalpy (surface and basal temperature) equations.

The sensitivity of total Antarctic ice volume to the initial ice thickness compares well with the calculated finite-difference-based value (Table 1, column 7), differing only by about 1 % at the fast-moving Slessor Ice Stream (Fig. 2a, location 3). Thickness sensitivities are relatively uniform and positive across the ice sheet over the 100-year simulation, except for a few outlet glaciers. A positive adjoint value in ice thickness indicates that a positive perturbation in ice thickness leads to a positive change in total volume, and vice versa. The Antarctic Ice Sheet is shown to have almost entirely positive adjoint values, as shown in Fig. 2a, except for a few marginal outlet glaciers. These few outlet glaciers that display negative ice thickness adjoint sensitivities contrast with other areas of ice discharge, notably the large floating ice shelves, which do not show any negative adjoint values. Figure 3a shows that the order of magnitude of this field of adjoint values is between 10⁸ and 10⁹ m², except for several very sensitive outlet glaciers, including Thwaites and Pine Island Glacier, glaciers in Marie Byrd, Oates, and Wilkes Land, and Byrd Glacier.

The pattern of the January (austral summer) precipitation adjoint values largely mirrors that of the ice thickness, with several distinctions. The order of magnitude is much larger, ranging instead between 10¹⁵ and 10¹⁷ m² (Fig. 3b). Outlet glaciers in Marie Byrd, Oates, Wilkes, and Queen Mary Land exhibit weaker sensitivities compared to the average Antarctic-wide summer precipitation sensitivities. Portions of floating ice shelves from Queen Maud eastward all the way to Queen Mary Land show the highest sensitivities overall to summer precipitation, while the larger ice shelves exhibit some of the lowest sensitivity to summer precipitation across the entire continent, almost an order of magnitude lower than the floating ice fringing the coast between Queen Maud and Queen Mary Land, from 10¹⁷ to 10¹⁶ m². Similar to the sensitivities to ice thickness, precipitation sensitivities are almost entirely positive, and the very lowest sensitivities are largely at the ice fronts (Fig. 3b). Table 1 shows less agreement in the January precipitation field calculated in the middle of the Ross Ice Shelf at point 2 (Fig. 2a), approximately a 12 % difference.

Sensitivities to surface and basal temperature (Fig. 2c and d) show much more structure than those of ice thickness and precipitation, with extended regions of positive and negative sensitivities. The sensitivities of total ice sheet volume to surface and basal temperature are largely negative, with the most negative values at the margins of the ice sheet and approaching zero toward the interior. The orders of magnitude of the surface and basal temperature sensitivities (Fig. 3c and d) are comparable to each other, with maximum values of approximately 10¹⁰ m^3 ∘C⁻¹.

Over the 100-year simulation, high sensitivities to surface and basal temperature at the margins extend inward toward the middle of the ice sheet following glacier drainage basins (Fig. 3c and d). There are two distinct differences between the sensitivities to surface and basal temperature seen in the adjoint fields. First, the highest sensitivities to basal temperature are higher than the highest sensitivities to surface temperature, indicating that the total ice sheet volume is in general more sensitive to changes in the applied basal temperature of the ice rather than at the surface in SICOPOLIS. Second, Fig. 3d shows that the location of those most sensitive areas to changes in basal temperature are at the grounding lines of ice shelves and glaciers, while the most sensitive areas to changes in surface temperature are all across the surface of the floating portions of ice, with the sensitivity increasing (becoming more negative) toward the ice fronts. The adjoint values of surface and basal temperature compare well with finite-difference-based sensitivities, differing by about 4 % and 8 %, respectively.

Table 1 also shows the results of a finite volume change calculation performed for the tile shown in Fig. 2a, region 4. Each control variable shown in column (1) is perturbed, as in the single point location finite differences, by ±5 % of the initial value of that field.

Columns (5) and (6) in Table 1 are computed in the following ways: $\mathrm{\Delta }{V}_{\mathrm{adj}}={\int }_{{\mathrm{\Omega }}_{\mathrm{4}}}\mathrm{\nabla }J\cdot \widehat{\mathit{n}}\mathrm{d}A$ is the area-integrated (or area-summed for discretized domains) volume change using the adjoint-derived gradient information, where Ω₄ is the area of the tile in Fig. 2a corresponding to the domain of control variable x. $\mathrm{\Delta }{V}_{\mathrm{fd}}=\frac{V(X+\mathit{\delta }X)-V(X-\mathit{\delta }X)}{\mathrm{2}}$ is the finite-difference-derived volume change, where X+δX is taken over the entire tile in region 4. In creating a sub-domain of Antarctica over which to calculate these finite volume changes, we selected an area of uniform sign. Areas with a great deal of sign variation might be more difficult to interpret since the adjoint values would tend to cancel each other.

The utility of this comparison is to convert the sensitivities into meaningful quantities that can be compared against each other to assess, for example, which control variable impacts the cost function the most given a perturbation of expected magnitude, in addition to providing another metric by which we may measure the adequacy of the adjoint model of SICOPOLIS. The percent difference between ΔV_adj and ΔV_fd over the 100-year time integration over tile 4 is largely higher than the point-wise measurements, ranging as high as 57 % for volume changes due to basal temperature perturbations uniformly in tile 4. Summer precipitation compares well, however, with a 6 % difference. Perhaps more interesting, the calculations in columns (5) and (6) suggest that, overall, summer precipitation has the largest impact on total ice sheet volume, with an approximate volume change of 10¹⁸ m³, compared against initial ice sheet thickness as well as surface and basal temperature perturbations, which at the lowest resulted in a volume loss of 10³ or 10⁴ m³ over 100 years. This result helps to explain the largest relative differences between adjoint and finite-difference sensitivities. These are large where the sensitivities are very small compared to the QoI, strongly suggesting that numerical noise plays an important role in either of these sensitivity calculations.

Lastly, the adjoint model of SICOPOLIS runs serially and completed 100 years of model runtime in 20, 75, and 600 min of wall clock time on a Linux box (Intel Xeon CPU E5-2650 at 2.00 GHz) for resolutions of 64, 40, and 20 km. The results shown in Figs. 2 and 3 are for 20 km resolution. To facilitate the adjoint computation SICOPOLIS-AD uses checkpointing, which is discussed in Appendix B.