1 Introduction↩︎

The Earth’s ice sheets, both present and past, are treated by the majority of ice sheet models as continuum materials that undergo viscous deformation described by the Stokes equations. Numerical models typically solve this set of momentum-balance equations, or some lower-dimensional approximation, for ice flow velocity in order to determine how ice is redistributed over time. In addition, many contemporary models contain machinery to compute gradients of model outputs with respect to certain input parameters - driven largely by a need to infer unknown basal or rheological properties of the ice using satellite-derived observations, for example, for the purpose of model initialisation [1]. For most models, this kind of inference is done considering the elliptic momentum balance problem in isolation (the case we focus on here) e.g. [2], [3]–[9], though some can accommodate time dependency as well [10]–[12].
There are two main methods for computing these gradients. The first is to rely on algorithmic differentiation (AD), should the model be built on a framework that allows it. This technology, used widely in the modelling of Earth system dynamics since the 1990s e.g. [13], [14]–[17], in theory, facilitates the computation of derivatives of complex functions expressed in code with little developer overhead. Programmatic operations are accompanied by rules defining the action of derivatives of those operations on tangents (known as forward mode AD) or cotangents (known as reverse or adjoint mode AD) to the parameter space. These rules can be derived before or at runtime, and there are a wide variety of techniques for doing so [18]. Various widely-used ice sheet models include some variant of this technology [10], [19], [20].
The alternative to AD, in the computation of first-order derivatives, is to use PDE-level adjoint methods. Here, gradients are found by solving a linear PDE for a set of adjoint variables which keep track of perturbations to the model outputs that maintain the validity of the model equations under input perturbations [21]. To remove some of the algebraic complexity, the first-order adjoint methods implemented in many models assume a linear rheology e.g. [2], [3], [8], [9], resulting in a self-adjoint set of model equations and approximate gradients. Of course, unless the approximate gradients are (roughly) parallel to the actual gradients, optimisation algorithms that depend on them will likely converge to different solutions. The use of the self-adjoint approximation has been found to produce inaccurate or low resolution solutions for scalar basal slipperiness parameters (widely used in ice sheet modelling) in various experiments involving synthetic data [22], [23]. However, the difference in performance when applied to real settings is thought to be much smaller [22], and has been shown, for some pan-Antarctic inverse problems, to be virtually zero [24]. This latter viewpoint has come to dominate over the last decade, and the widespread use of the linear-viscosity approximation in constructing first-order adjoints reflects the rarity of major errors arising from it.
Beyond the computation of gradients, there is growing interest in second-order derivative information among ice sheet modellers. For example, Hessian-vector products enable certain powerful optimisation methods (e.g. Newton’s method), and are essential for many Bayesian frameworks for parametric uncertainty quantification [2], [25]–[29]. To date, codes involving second-order derivatives in the ice sheet modelling community have been based on finite-element libraries (like FEnICs/Firedrake) employing a particular type of AD applicable to that domain [12], [20].
Models lacking AD can, once again, use PDE-level adjoint methods to compute these second-order derivatives - at higher accuracy and efficiency than methods based, at least in part, on finite differences (e.g. [2]). Extending the logic of the first-order method described above, second-order derivatives can be computed by constructing two further adjoint variables to keep track of: 1) second-order changes to the state variable and 2) first-order changes to the first-order adjoint variable (see, e.g. [30], [31]). This incurs two additional linear solves. In principle, this can be repeated ad infinitum, computing \(n(n+1)/2\) adjoint variables, using \(n(n+1)/2\) linear solves, in pursuit of the \(n^{\rm th}\)-order derivative. The barrier to implementation is not so much the computational overhead, but the algebraic one in deriving the equations each adjoint variable must obey. However, pockets of the ice sheet modelling community have utilised complete second-order adjoint methods, formulated using variational methods, to good effect: in performing inverse problems with finite element full-Stokes solvers, demonstrating the efficacy of the method in improving the efficiency of inverse solvers and its utility in uncertainty quantification [32]–[34].
Here, we derive a PDE-level second-order adjoint formulation for the shallow-stream approximation (SSA) to the momentum balance, using a linear viscosity approximation to mirror common practice in first-order adjoint methods. We label this the second-order self-adjoint (SOSA) method. The resulting formulation can, therefore, be straightforwardly integrated into a variety of existing SSA (or SSA-like) models and used to compute approximate Hessian-vector products. We implement this model in a finite volume code written in the numerical computation library JAX [35]. The model is similar to a uniform-mesh version of the BISICLES ice sheet model [36], but fully differentiable, allowing us to compare the SOSA Hessian with an exact Hessian obtained via AD. We qualitatively compare Hessian-vector products (HVPs) in different cases, and perform spectral analysis on each Hessian. We find that, while HVPs from the two methods appear qualitatively similar, the invariant subspaces spanned by their leading eigenvectors begin to diverge almost immediately. However, this divergence ceases after around 50 modes, suggesting a diminishing importance in the neglected nonlinearity. We judge the full Hessian to be preferable where possible, though there may be situations in which SOSA can provide an adequate approximation, with fidelity that depends on the rank of the approximation.

2 The adjoint formulation↩︎

First, we briefly outline the first- and second-order adjoint models for a generic numerical model, before specifying the equations for the SSA formulation of the momentum balance under the self-adjoint approximation. For a full derivation, see 8 and 9.
Let \(U\) and \(Q\) be Hilbert spaces with the standard L2 inner product \(\langle\cdot,\cdot\rangle\). Let us assume a model \(\mathcal{M}\) with input parameters \(q\in Q\) and outputs \(\mathbf{u}\in U\), such that \(\mathcal{M}: q\mapsto \mathbf{u}\), governed by the system of partial differential equations \(G(\mathbf{u}, q)=0\), and a functional of interest \(\mathcal{J}:U\times Q\to \mathbb{R}\). This could be, for example, the total grounding line flux over some portion of an ice sheet or a cost function comparing modelled and observed ice velocity. The quantities of interest are the gradient \(d_q\mathcal{J}:T_Q\to \mathbb{R}\) and the Hessian \(d^2_q\mathcal{J}= \mathcal{H}_q:T_Q\times T_Q\to \mathbb{R}\), where \(T_Q\) is the tangent space to \(Q\) at the current parameter value.
We write the action of the gradient on a perturbation \(\delta q\) as \(d_q\mathcal{J}[\delta q]=\langle d_q\mathcal{J}, \delta q\rangle\) and the action of the Hessian as \(d^2_q\mathcal{J}[\delta q_1, \delta q_2] = \langle d^2\mathcal{J}[\delta q_1], \delta q_2\rangle\). The Hessian-vector-product \(d^2_q\mathcal{J}[\delta q_1]\) is itself a gradient operator. The notation here relies on the fact that \(T_q\cong Q\), \(T_\mathbf{u}\cong U\) and the Reisz representation theorem to associate gradients with elements of \(U\) and \(Q\).

2.1 General first- and second-order adjoint models↩︎

The total variation of the functional \(\mathcal{J}\) under a perturbation \(\delta q\) has, according to the chain rule, two components: \[\label{dJ} d_q\mathcal{J}[\delta q] = (\partial_q\mathcal{J}+ \partial_\mathbf{u}\mathcal{J}~d_q\mathbf{u})[\delta q].\tag{1}\] So we write \(d_q\mathcal{J}= \partial_q\mathcal{J}+ \partial_\mathbf{u}\mathcal{J}~d_q\mathbf{u}\). The adjoint method is a handy trick for calculating this without explicitly computing the Jacobian \(d_q\mathbf{u}\), which would be prohibitive. Under a small perturbation \(\delta q\) that preserves the properties of the original parameter field, the solver will converge just as well, but to a slightly different \(\mathbf{u}\). Hence: \[\begin{align} d_q\mathcal{J}[\delta q] &= d_q\mathcal{J}[\delta q] + d_q G[\delta q][\boldsymbol{\lambda}]~~~~~~ \forall~\boldsymbol{\lambda}\in U \end{align}\] because \(d_q G[\delta q][\boldsymbol{\lambda}] = 0~~~ \forall~\boldsymbol{\lambda}\in U\) .
The trick is to choose the adjoint variable \(\boldsymbol{\lambda}\) such that this formula contains no reference to the troublesome Jacobian \(d_q\mathbf{u}\). Following it through (8), one finds that if \(\boldsymbol{\lambda}\) solves the following system of PDEs: \[\label{foa95system} \partial_\mathbf{u}\mathcal{J}+ (\partial_\mathbf{u}G)^\ast[\boldsymbol{\lambda}] = 0,\tag{2}\] then: \[\label{foa95grad} d_q\mathcal{J}= \partial_q\mathcal{J}+ (\partial_q G)^\ast [\boldsymbol{\lambda}].\tag{3}\] The equations 2 define the First Order Adjoint (FOA) system; so-called because \((\partial_{\mathbf{u}}G)^\ast\) is the adjoint of the tangent linear operator. The parameter \(\boldsymbol{\lambda}\) encodes the changes to \(\mathbf{u}\) required to satisfy the model equations under perturbations of the input parameters.
Many applications require only that we can compute Hessian-vector products of the form \(d^2_q\mathcal{J}[\delta q]\), rather than an explicit representation of the full Hessian. The second order adjoint method allows us to do this efficiently by applying the same trick as the first order method, this time defining an additional two adjoint variables \(\boldsymbol{\mu},\boldsymbol{\beta}\in U\). These, respectively, keep track of second order changes to \(\mathbf{u}\) under perturbations required to maintain the validity of the model equations, and first order changes to the adjoint variable \(\boldsymbol{\lambda}\) required to maintain the validity of the FOA equations. The derivation is presented in 9.
Let us define: \[G_{\boldsymbol{\lambda}}= \langle\boldsymbol{\lambda}, G\rangle\] Then the Hessian-vector product can be written: \[\label{generic95hvp} d^2_q\mathcal{J}[\delta q] = (\partial_{q}^2\mathcal{J}+ \partial_{q}^2G_{\boldsymbol{\lambda}})[\delta q] + (\partial_{q}G)^\ast[\boldsymbol{\beta}] + (\partial_{\mathbf{u}}\partial_{q}\mathcal{J}+ (\partial_{q}G_{\boldsymbol{\lambda}})^\ast)[\boldsymbol{\mu}]\tag{4}\] where \(\boldsymbol{\mu}\) and \(\boldsymbol{\beta}\) solve a set of second-order adjoint (SOA) equations of the form: \[\begin{align} &\partial_{q}\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\delta q] + \partial_{\mathbf{u}}\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\boldsymbol{\mu}] = 0\tag{5}\\ &(\partial_{q}\partial_{\mathbf{u}}\mathcal{J}+ \partial_{q}\partial_{\mathbf{u}}G_{\boldsymbol{\lambda}})[\delta q] + (\partial_{\mathbf{u}}G)^\ast[\boldsymbol{\beta}] + (\partial_{\mathbf{u}}^2\mathcal{J}+ \partial_{q}^2G_{\boldsymbol{\lambda}})[\boldsymbol{\mu}] = 0.\tag{6} \end{align}\]

2.2 Application to the shallow stream equations↩︎

The FOA and SOA systems given above are generic to any model solving some system of equations expressible as \(G(\mathbf{u}, q)=0\). In this study, we look specifically at those which solve the shallow stream approximation (SSA) formulation of the momentum balance governing an ice stream [37]. We shall assume without loss of generality that the stress at the base of the ice is a linear function of basal ice velocity. We focus on gradients with respect to the stiffness of the ice, and this is trivially extended to include the slipperiness.
Defining the scalar field \(\varphi(q)=\phi_0 e^{q}\) and the resistive stress tensor field \(H(\mathbf{u}) = h\bar{\mu}\left[\grad\mathbf{u}+(\grad\mathbf{u})^\top+2(\grad\cdot\mathbf{u})\mathcal{I}\right]\), the SSA equations read:

\[\label{eq:ssa} G(q, \mathbf{u}) = \grad\cdot[\varphi(q)H(\mathbf{u})] - C\mathbf{u}- \rho_i gh\grad s = 0,\tag{7}\]

where \(h\) is the thickness of the ice, \(\mathbf{u}\) is the two-dimensional ice velocity, \(\bar{\mu}\) is the vertically-averaged effective viscosity, \(C\) is a basal sliding coefficient, \(\rho_i\) is the density of ice, \(g\) is the acceleration due to gravity, and \(s\) is the height of the surface of the ice above sea level.
For any \(\boldsymbol{\lambda}, G(\mathbf{u}, q) \in U\) and \(\delta q \in Q\), we can write: \[\label{obvs95identity} (\partial_{\mathbf{u},q} G)^\ast[\boldsymbol{\lambda}] = (\partial_{\mathbf{u},q} \langle \boldsymbol{\lambda}, G \rangle).\tag{8}\] Given this, we can write the FOA system as: \[\label{foa95system95ssa} \partial_\mathbf{u}\mathcal{J}+ \partial_\mathbf{u}G_{\boldsymbol{\lambda}}= 0,\tag{9}\] and: \[d_q\mathcal{J}= \partial_q\mathcal{J}+ \partial_q G_{\boldsymbol{\lambda}}\] where \(G_{\boldsymbol{\lambda}}=\langle \boldsymbol{\lambda},G(\mathbf{u},q) \rangle\), as above.

2.2.1 First order adjoint model↩︎

We consider cases in which the scalar and vector fields of interest and their spatial gradients vanish on the boundaries of the domain. This restricts \(U\) and \(Q\) from being generic Hilbert spaces. With this, there is a natural equivalence between taking the adjoint of an operator and performing integration by parts. It is easy to ensure this with the finite volume method. When using finite elements, the boundary conditions imposed on \(\mathbf{u}\) the edges of the domain must also be imposed on the adjoint variables.
With this, the first order adjoint system for the SSA equations is: \[\label{foa95ssa} \nabla\cdot(\varphi H(\boldsymbol{\lambda})) - C\boldsymbol{\lambda}+ \partial_{\mathbf{u}}\mathcal{J}= 0.\tag{10}\]

With \(\boldsymbol{\lambda}\) as the solution to this boundary value problem, the gradient can be written: \[\label{gradient95ssa} d_q\mathcal{J}= \partial_{q}\mathcal{J}- \nabla\boldsymbol{\lambda}:\varphi H(\mathbf{u}).\tag{11}\]

2.2.2 Second order adjoint model↩︎

There are 12 terms to consider in the SOA model and the Hessian-vector product defined in equations (4 6 ). However, once more making the assumption of linearity in the viscosity, their form becomes familiar again. The SOA system for the SSA equations is the following set of elliptic PDEs:

\[\begin{align} &\grad\cdot(\varphi\delta qH(\mathbf{u})) + \nabla\cdot(\varphi H(\boldsymbol{\mu})) - C\boldsymbol{\mu}= 0\tag{12}\\ &\grad\cdot(\varphi\delta qH(\boldsymbol{\lambda})) + \nabla\cdot(\varphi H(\boldsymbol{\beta})) - C\boldsymbol{\beta}+ (\partial_{\mathbf{u}}^2\mathcal{J}, \boldsymbol{\mu})= 0\tag{13}. \end{align}\] We have made the assumption here that the functional of interest is separable in functions involving \(\mathbf{u}\) and \(q\) so that \(\partial_{q}\partial_{\mathbf{u}}\mathcal{J}=0\).
With \(\boldsymbol{\mu}\) and \(\boldsymbol{\beta}\) found by solving the above, the Hessian vector product is:

\[\begin{align} \label{hvp95final} (\mathcal{H},\delta q) = (\partial^2_q\mathcal{J},\delta q) - \nabla\boldsymbol{\lambda}:(\varphi\delta q H(\mathbf{u})) -\nabla\boldsymbol{\lambda}:\varphi H(\boldsymbol{\mu}) -\nabla\boldsymbol{\beta}:\varphi H(\mathbf{u}). \end{align}\tag{14}\]

This shows that, under the self-adjoint approximation, the first- and second-order systems should be solvable by any code built for the SSA equations. We see the same terms repeatedly appearing: the gradient of something that looks like the stress tensor, the linear sliding term and a right-hand side that depends on the form of our functional.

3 The algorithmic differentiation formulation↩︎

Starting with equation 1 , algorithmic differentiation provides an alternative set of methods for dealing with the term \(\partial_{\mathbf{u}}\mathcal{J}d_q\mathbf{u}\). A simple method, implemented, for example, by neural networks, is to unroll the set of elementary operations \(\{f_i\}\) that make up the mapping between \(q\) and \(\mathbf{u}\). Imagine \(\mathbf{u}= f_n\circ f_{n-1}\circ\cdots\circ f_0(q)\), then: \[(\partial_{\mathbf{u}}\mathcal{J}~d_q\mathbf{u})|_q = (\partial_{\mathbf{u}}\mathcal{J}|_\mathbf{u})(\partial_{f_{n}}\mathbf{u}|_{f_{n}})(\partial_{f_{n-1}}f_{n}|_{f_{n-1}})\cdots(\partial_{q}f_0|_{q}),\] (though, in all likelihood, there will be a more complicated computational graph than this structure shows, probably consisting of many branches). In adjoint mode AD (also known as ‘back-propagation’) the series of operations is performed left-to-right. This is significantly more computationally efficient than performing the operations right-to-left when \(f_i\) are high-dimensional and the final operation is a functional. However, in a numerical model, the number of operations is huge, and stacking them in this way is not possible.
Instead of this, we find a more computationally tractable approach is to implement equations 2 and 3 but use AD to define the action of \((\partial_{\mathbf{u}}G)^\ast\) and \((\partial_{q}G)^\ast\) on adjoint variables \(\boldsymbol{\lambda}\) fully - including the nonlinear dependence of ice viscosity on the speed \(\mathbf{u}\). The function for computing these gradients can then, itself, be differentiated using AD. We use adjoint mode AD for this task, meaning that functions must be specified for the propagation of adjoint variables through a linear solver, but this is again a simple application of a first-order adjoint method.

4 Implementation↩︎

We make use of a numerical computation library for Python called JAX [35] to implement the SSA model (eq. 7 ) and its various adjoints (equations 10 14 ). With this, we solve the problems on a uniform mesh using the finite volume method. For the forward problem, we use Newton’s method, and utilise JAX’s AD functionality to compute the Jacobian of the residual equations during each iteration. We use PETSc’s python interface [38] to solve linear algebra problems directly using MUMPS [39]. To compute the spectral decompositions of the Hessians, we used the Lanczos method using SciPy’s sparse.linalg.eigsh function [40]. Throughout, we use a Glen’s flow law exponent of \(n=3\).

4.1 Test domains↩︎

To compare the Hessians produced using the AD and SOSA methods, we constructed two synthetic domains: one representing an ice shelf and the other an ice stream (figure 1). The first (labelled “Ice Shelf” - figure 1a) consists of a slab of ice of uniform \(500~{\rm m}\) thickness, and a slipperiness coefficient \(C\) defined to be zero under most of the ice, and high (\(C=10^4\)) in a narrow band around three of the sides. The calving front is situated two cells in from the right of the domain. We applied reflection boundary conditions to the velocity at each boundary, though the high \(C\)-value at three of the boundaries, and the lack of ice at the other, means that the velocity and its gradients are zero there.
The other domain (labelled “Twisty Stream” - figure 1b) is borrowed from [36]. It represents an infinitely long, uniform-thickness (\(1~{\rm km}\)) ice stream following the line of a snaking valley. The slipperiness coefficient is defined as \[C = 10^3\times (1 + \epsilon + \sin{\left(\frac{\pi}{2} + 2\pi\frac{y}{R} + m\sin{(2\pi \frac{x}{R})}\right)}\] where \(\epsilon=5\times10^{-3}\), \(R=180~{\rm km}\), and \(m=1/4\). The bed is sloping to the right with uniform slope of \(0.5^\circ\). To this, we applied periodic boundary conditions on the the left- and right-hand boundaries, and reflection boundary conditions at the top and bottom.
The speed solutions for these domains are shown in figure 2a and 2e, respectively.

5 Results↩︎

5.1 Direct comparisons↩︎

First, we show a visual comparison between example Hessian-vector products found using the SOSA and AD methods. We considered the functional \(\mathcal{J}(\mathbf{u}) = \int_\Omega\sqrt{\mathbf{u}\cdot\mathbf{u}}~{\rm d}\Omega\). For both the Ice Shelf and Twisty Stream, we computed the ice speed, the gradient of \(\mathcal{J}\) with respect to the parameter \(q\) (defined in equation 7 ), then the Hessian-vector product \(d^2_q\mathcal{J}[\delta q]\) where the perturbation \(\delta q= \epsilon d_q\mathcal{J}\) is in the gradient direction. Results are shown in Fig. 2.
In the ice shelf case, the gradient reflects the common observation that softening in the shear margins and in a central compressive band [41] have largest impact on flow speed. The HVP in this direction, for both the AD and SOSA cases, display a similar structure. Both indicate a decreasing sensitivity to softening as it is applied upstream of the compressive arch, but an increasing sensitivity to it downstream. In the AD case, the magnitudes of the HVP are larger by up to an order of magnitude.
The Twisty Stream experiment shows that the ice speed, as expected, is most sensitive to softening in the shear margins of the ice stream. In both the AD and SOSA cases, the sensitivity to softening increases as it is applied at the outside edge of the current shear margins (corresponding to a widening of the ice stream and an increase in integrated flow speed), and a decreased sensitivity to softening at the inside edge (corresponding to a narrowing of the ice stream). Once more, this effect is amplified when the full adjoint is found using AD than in the SOSA case.

5.2 The Hessian decomposition↩︎

To compare the methods more concretely, we reduce our attention to the Twisty Stream example and consider the spectral decomposition of the Hessians computed via AD and SOSA. The domain is \(180\times 180\)-cells in size, so the Hessian has over 32 000 eigenvectors. We consider only the first 500 to get an idea of how the Hessians differ in their most consquential modes.
Given that we have in mind uses of the Hessian in the context of glaciological inverse problems, we look at a functional representing something like the misfit between modelled and observed ice speed: \[\label{cst95fct95ip} \mathcal{J}= \int_\Omega \sqrt{(\mathbf{u}-\mathbf{u}_{\rm sltn})\cdot(\mathbf{u}-\mathbf{u}_{\rm sltn})}~{\rm d}\Omega\tag{15}\]

where \(\mathbf{u}_{\rm sltn}\) is the ice velocity solution for unperturbed \(q\). This has the feature of having zero gradient at \(\mathbf{u}= \mathbf{u}_{\rm sltn}\).

Fig 3 shows the first 10 eigenvectors produced using each method for this functional. Many similar patterns appear in the two sets, but their order is different, and often differ by the sign of the eigenvalue (e.g. the principal mode in each). For example, the pair of eigenvectors 4 and 5 of the AD set look similar to the 6-7 pair of SOSA eigenvectors. This pattern continues throughout the first 500 eigenvectors, with similar structures appearing in both sets. There is also a difference in the along-flow/across-flow periodicity in some early eigenvectors. For example, eigenvectors 4 and 5 of the SOSA-derived Hessian show half-wavelength in the across-flow direction and 2 wavelengths in the along-flow, and no such pattern appears in the first 500 eigenvectors of the AD-derived Hessian.
The eigenvalues for each Hessian reduce rapidly in magnitude for the first few modes, though go onto display sub-exponential decay (Fig. 4). This behaviour has been seen in the highest order eigenvalues for similar problems before, and has been shown to give way to exponential decay as the order gets large [28]. The SOSA eigenvalues are smaller than the AD eigenvalues, due to the missing terms associated with derivatives of the viscosity. For the first few eigenvalues, the ratio is around a factor of 3 (Fig. 4 - blue line). This relates to the fact that the operator \(\partial_{\mathbf{u}}G\) is roughly a factor of \(n=3\) (our choice for Glen’s exponent) smaller when derivatives of the viscosity are included. This factor of \(n\) propagates unchanged into the first term of equation 14 and into the eigenvalues for the eigenvectors in which this term dominates. After the first 5 modes, the difference drifts away from a constant factor, though remains largely between 3 and 4, though will drop below 3 eventually.

The similarity of the sets of eigenvectors under visual inspection means that greater analysis is needed to determine the utility of the SOSA approximation. For most applications in optimisation and uncertainty quantification, an approximation of the Hessian to some rank much smaller than the full dimension is necessary. For example, in an idealised optimisation procedure, one can imagine the search directions being described by the ordered eigenvectors of the Hessian of the cost function, and we don’t need to proceed through all of them before being satisfied the problem is solved well-enough. So, we need to determine the similarity of the AD and SOSA Hessians as the rank of the approximation increases. Concretely, we ask at what point the subspaces spanned by the first k eigenvectors of each Hessian start to diverge.
To determine this, it is necessary to first consider whether the individual subspaces are themselves self-consistent, so that any divergence between the two reflects the true geometry rather than, for example, numerical error. One might worry about this, for example, because the methods include a number of linear algebra solves per HVP with poorly conditioned matrices. For each eigenvector/eigenvalue pair \((\mathbf{\hat{e}}_i, c_i)\), we computed an eigenvalue residual: \[\label{eigres} {\rm eigres_i} = \frac{||d_q^2\mathcal{J}[\mathbf{\hat{e}}_i] - c_i\mathbf{\hat{e}}_i||_1}{||c_i\mathbf{\hat{e}}_i||_1}.\tag{16}\] Figure 5b shows these remain small throughout the first 500 pairs computed for each method. Additionally, we tested the orthonormality of the subspaces by computing, for each subspace, the metric: \[\label{orthres} {\rm orthres_k} = ||S^\top_k S_k - \mathcal{I}_k||_1,\tag{17}\] where \(S_k = [\mathbf{\hat{e}}_1, \dots,\mathbf{\hat{e}}_k]\) and \(\mathcal{I}_k\) is the \(k\times k\) identity matrix. Figure 5a shows that this too remains small up to \(k=500\). The conclusion is that each method indeed produces a clean set of orthonormal eigenvectors.
An intuitive notion of the similarity of subspaces constructed from the leading modes of the AD and SOSA Hessians can be determined by considering the set of principal (or canonical) angles between them. Let \(A\) be the matrix with columns defined by the first \(k\) eigenvectors found using the SOSA method and \(B\) be the equivalent matrix for the AD method. The singular value decomposition of the matrix \(A^\top B\) defines the cosines of the principal angles between the subspaces defined by the columns of \(A\) and \(B\). When the angles are all close to zero, the projection of any vector in one subspace onto the other will be almost parallel to the first. When the largest of these angles is \(\pi/2\), the subspaces are orthogonal, though the projection distance still might be small. As such, the proportion of angles significantly different from 0 gives us an indication of the similarity of the Hessians at different rank. Figure 5c shows the ordered principal angles between subspaces constructed using the first \(k=1\,{\rm to}\,500\) eigenvectors of the two Hessians.

Figure 5c shows that the subspaces spanned by the eigenvectors of the SOSA and AD Hessians start out close to parallel but become persistently orthogonal after the 33rd mode. Up to the first 4 modes, the principal angles are all small (\(<\pi/8\)), however, after this, the proportion of larger angles grows until around the 50th mode (5c - inset). From here, the proportion of angles very close to 0 (white region) starts to dominate again. The number of angles near \(\pi/2\) (black region) grows as modes are added, but it does not increase as a proportion of the full set of angles. Hence, though there is a persistent inaccessible region to the SOSA Hessian, increasing its rank appears to always improve the approximation. This general pattern implies that there is a set of ranks for which neglecting nonlinear terms makes the SOSA Hessian most different to the full Hessian: between \(4\) and \(\sim100\). This analysis is revealing about why gradient methods that assume self-adjointness appear to work well in practice, and has consequences for the use of the SOSA method in place of the full Hessian. We discuss these points below.

6 Discussion↩︎

6.1 Optimisation↩︎

As mentioned above, the context we have in mind for using the Hessian is the solution of inverse problems, for example for the purposes of model initialisation or diagnostic modelling, and quantifying uncertainty in the solutions. Commonly, solutions are found using gradient-based methods using only first order information. Despite not explicitly computing second-order derivatives, the geometry of the cost-landscape (e.g. 15 ) encoded in the Hessian still governs how those methods behave. We imagine, as an approximation, that reasonably close to the solution, steepest-descent-type methods (e.g. conjugate gradients) explore in turn the maximum-curvature directions corresponding to the eigenvectors of the Hessian. During an iterative optimisation algorithm using an approximate gradient, the solution will continue to improve only while the search direction remains non-orthogonal to the true gradient. The plot of principal angles between the AD- and SOSA-derived Hessians (Fig. 5c) indicates that these directions are similar in the first few modes, then diverge slightly for the next few tens of modes, before divergence stops around mode 100. This gives us some insight into why gradient-based methods applying the self-adjoint approximation perform well in practice - there is no obvious point at which adding modes corresponding to inexact gradients stops improving the solution. Because they are expensive for real problems, optimisation algorithms are, in practice, run only for a few tens of iterations. However, Fig. 5 suggests improvement would continue on from this.
One potential use of the SOSA Hessian is in higher-order Newton-like methods, where we could use it to approximate the true Hessian at some low rank. The same considerations apply here as in the first-order case - in that the invariant subspaces spanned by the eigenvectors of each Hessian are similar enough to imagine this being effective. However, underestimation of curvature in each mode of the SOSA Hessian now becomes important. Depending on the rank to which the approximation is desired, a scaling factor could be applied to the eigenvalues (Fig. 4). Given the difference in many of the leading-order curvature modes, it is unclear how well even a scaled SOSA Hessian would perform. Hence, where feasible, the full Hessian is preferable.

6.2 Uncertainty Quantification↩︎

Finally, in Bayesian parametric uncertainty quantification, the inverse Hessian is used to calculate the covariance of the posterior distribution for the parameter of interest [12], [26], [33]. The accuracy of the posterior covariance therefore places strict requirements on the fidelity of the Hessian. Given that the eigenvalues are different in the two cases, use of the SOSA Hessian would result in a systematic underestimation of uncertainty in each direction. To some degree, it seems that this could be compensated-for in the leading modes by applying a constant scaling to the inverse Hessian, though this would be quite inexact. Existing studies have shown the number of eigenvectors required for faithful construction of the posterior covariance numbers in the thousands, for problems not much larger than those considered here [28]. Even though the proportion of large angles between the SOSA and full Hessian might be small at this order, the differences at low order, and differences in eigenvalues, might well yield inaccurate results. However, there exist other methods of uncertainty quantification, such as perturbation ensembles, for which the broad alignment of the SOSA and full Hessian might suffice to produce similar results.

7 Conclusion↩︎

The increasing prevalence of algorithmic differentiation tools has led to the adoption of methods in ice sheet model initialisation that make use of second-order derivative information. The second-order adjoint model derived here, which neglects terms coming from the nonlinear dependency of ice viscosity on velocity, is an easily implemented alternative for the set of models not built on frameworks which support AD. The Hessian-vector products produced by this SOSA method, implemented in a JAX-based finite volume code, look, in the two cases we consider, to be visually similar, as do the first 500 eigenvectors of the Hessian in each case. Spectral analysis reveals that there is divergence between the subspaces spanned by eigenvectors of the full and approximate Hessians. This divergence is small enough to explain the success of first-order optimisation methods using the self-adjoint approximation. However, coupled with the difference in eigenvalues, it suggests use of the full Hessian, derived via AD, ought to be used where a high-fidelity approximation is sought at any rank higher than 3 or 4.

Acknowledgements TSS and SLC are supported by the Advanced Research and Invention Agency grant no. SCOP-PR01-P010.

8 First-order Adjoint Formulation↩︎

We consider the sensitivities of a functional \(\mathcal{J}(\mathbf{u},q)\in \mathbb{R}\) with respect to changes in the control parameter \(q\), for a model that solves a set of partial differential equations \(G(\mathbf{u},q)=0\). We take \(q\in Q\) and \(\mathbf{u}\in U\) where \(Q\) and \(U\) are both Hilbert spaces, so we use inner products throughout this derivation, and apply the Reisz representation theorem wherever convenient to associate gradients of functionals and Hessian-vector products with elements of those different spaces.
The ultimate quantity of interest in calculating the sensitivity of \(\mathcal{J}\) changes in \(\phi\), is the Gâteaux derivative of \(\mathcal{J}\) in the direction \(\delta q\): \[\label{gateaux95def} \delta\mathcal{J}(\mathbf{u}, q, \delta q) = \lim_{\varepsilon\rightarrow0}\left\{\frac{1}{\varepsilon}\left(\mathcal{J}(\mathbf{u}, q+\varepsilon\delta q) - \mathcal{J}(\mathbf{u}, q)\right)\right\}.\tag{18}\]

We define the gradient \(d_q\mathcal{J}\) as an operator, and as an element of \(Q\) by equating all of the following things: \[\label{eq:grad95direction95ip} \delta\mathcal{J}(\mathbf{u}, q, \delta q) = d_q\mathcal{J}[\delta q] = \langle~d_{q}\mathcal{J}, \delta q~\rangle = \int_\Omega d_q\mathcal{J}(\mathbf{x})\, \delta q(\mathbf{x})~d\Omega.\tag{19}\]

The gradient \(d_q\mathcal{J}\) can be decomposed via the chain rule, which operates in the same way in functional calculus as calculus:

\[d_q\mathcal{J}= \partial_{q}\mathcal{J}+ \partial_{\mathbf{u}}\mathcal{J}\,d_q\mathbf{u}.\]

The first-order adjoint method allows us to compute this without forming the Jacobian \(d_q\mathbf{u}\) explicitly. There are a few ways of formulating it, but it is typical in the ice sheet modelling literature to imagine forming the Lagrangian: \[\mathcal{L}= \mathcal{J}+ \langle \boldsymbol{\lambda}, G\rangle\] for some Lagrange multipliers \(\boldsymbol{\lambda}\in U\). These Lagrange multipliers are also known as the adjoint variables. The variation of this under perturbations of the control parameter \(q\) is simply the variation of \(\mathcal{J}\), given that \(d_q G[\delta q]=0\). (In other words, we still expect our solver to converge when \(q\) is perturbed a small amount.) So, let us consider the action of the gradient of \(\mathcal{L}\) on a perturbation \(\delta q\):

\[\label{lag95grad} d_q\mathcal{J}[\delta q] \equiv d_q\mathcal{L}[\delta q] = \big(\partial_{q}\mathcal{J}+ (\partial_{\mathbf{u}}\mathcal{J}\,d_q\mathbf{u}) + \partial_{q}\langle\boldsymbol{\lambda},G\rangle + (\partial_{\mathbf{u}}\langle\boldsymbol{\lambda},G\rangle\,d_q\mathbf{u})\big)[\delta q].\tag{20}\]

From the definition of the Gateaux derivative 18 , it can be seen that:

\[\label{obvs95identity95back} (\partial_{\mathbf{u},q} \langle \boldsymbol{\lambda}, G \rangle)[\delta q] = \langle\boldsymbol{\lambda},\partial_{\mathbf{u},q} G[\delta q]\rangle,\tag{21}\] In which case, 20 gives us:

\[d_q\mathcal{J}= \partial_{q}\mathcal{J}+ (\partial_{q}G)^\ast[\boldsymbol{\lambda}] + (\partial_{\mathbf{u}}\mathcal{J}+ (\partial_{\mathbf{u}}G)^\ast[\boldsymbol{\lambda}])d_q\mathbf{u}\] where we have used the definition of the adjoint of a linear map.
Hence, if we can solve the first-order adjoint equation: \[\label{generic95foa95eq} \partial_{\mathbf{u}}\mathcal{J}+ (\partial_{\mathbf{u}}G)^\ast[\boldsymbol{\lambda}] = 0\tag{22}\] for \(\boldsymbol{\lambda}\), then: \[\label{grad95given95foa} d_q\mathcal{J}= \partial_{q}\mathcal{J}+ (\partial_{q}G)^\ast[\boldsymbol{\lambda}].\tag{23}\]

The question now becomes: what what are \((\partial_{\mathbf{u}}G)^\ast[\boldsymbol{\lambda}]\) and \((\partial_{q}G)^\ast[\boldsymbol{\lambda}]\) for the SSA equations?

8.1 The adjoint equation and gradient for the SSA↩︎

Once, again the shallow-stream approximation equations are: \[\grad\cdot[\varphi(q)H(\mathbf{u})] - C\mathbf{u}- \rho_i gh\grad s = 0.\] where \(\varphi(q) = \phi_0e^q\) and \(H(\mathbf{u})\) is the resistive stress tensor, subject to boundary conditions: \[\label{eq:adj95bvp952} \begin{align} \hat{\mathbf{n}}\cdot\mathbf{u}= 0,\quad \hat{\mathbf{t}}\cdot\grad\mathbf{u}\cdot\hat{\mathbf{n}} = 0 \end{align}\tag{24}\] (where \(\hat{\mathbf{n}}\) and \(\hat{\mathbf{t}}\) are normal and tangent vectors to the boundary respectively).
The linear self-adjoint approximation, used widely in the literature, can be simply stated as: \[\label{linsa} (\partial_{\mathbf{u}}G(\mathbf{u}))^\ast[\boldsymbol{\lambda}] = G(\boldsymbol{\lambda}).\tag{25}\]

Given the structure of \(G\), it is clear that making \(\bar{\mu}\) independent of \(\mathbf{u}\) makes \(G\) linear. In order to justify the self-adjointness, one must ensure that adjoint variables \(\boldsymbol{\lambda}\) also conform to the boundary conditions specified in 24 . The result of this is that all the boundary terms arising from integration by parts disappear when deriving 25 .
The \((\partial_{q}G)^\ast[\boldsymbol{\lambda}]\) term is most easily computed by rewriting \((\partial_{q}G)^\ast[\boldsymbol{\lambda}] = \partial_{q}\langle\boldsymbol{\lambda}, G\rangle[\delta q]\) and then direct application of the definition of the Gateaux derivative 18 . Briefly: \[\begin{align} (\partial_{q}G)^\ast[\boldsymbol{\lambda}] &= \partial_{q}\langle\boldsymbol{\lambda}, G\rangle[\delta q]\\ &= - \lim_{\varepsilon\rightarrow0}\left\{\frac{1}{\varepsilon}\left(\int\nabla\boldsymbol{\lambda}: ([\varphi(q+\varepsilon\delta q) -\varphi(q)]H)~d\Omega)\right)\right\}\\ &= -\nabla\boldsymbol{\lambda}: (\varphi H). \end{align}\] Where, the second line makes use of the following statement, which is just generalisation of Green’s first identity: For an n-dimensional vector field \(\mathbf{v}\) and a rank-2 tensor field \(T\) defined over the domain \(\Omega\): \[\int_{\Omega}\mathbf{v}\cdot(\nabla\cdot T)~d\Omega = \int_{\partial\Omega}(T~\cdot\hat{\mathbf{n}})\mathbf{v}~d\Sigma - \int_{\Omega}T\colon(\nabla\mathbf{v})~d\Omega,\] where \(\hat{\mathbf{n}}\) is the unit normal to the boundary \(\partial\Omega\).
With this, we find: \[\label{eq:gradient} d_q\mathcal{J}\approx \partial_{q}\mathcal{J}- \grad\boldsymbol{\lambda}:(\varphi(q)H(\mathbf{u})).\tag{26}\]

9 Second-order Adjoint Formulation↩︎

9.1 General formulation↩︎

We start by restating the adjoint equations we found in the derivation of the first-order sensitivity: \[\label{adjoint95eq952} \partial_{\mathbf{u}}\mathcal{J}+ \partial_{\mathbf{u}} G_{\boldsymbol{\lambda}} = 0.\tag{27}\] These are still accompanied by the constraint that \(\mathbf{u}\) solves the stress-balance equations: \[\label{constraint952} G=0.\tag{28}\] The first step in calculating second-order derivatives is to solve the first order equations for the Lagrange multipliers \(\boldsymbol{\lambda}\). This allows us to write the functional gradient of the cost function as: \[d_q\mathcal{J}= \partial_{q}\mathcal{J}+\partial_{q} G_{\boldsymbol{\lambda}}.\]
Differentiating this and applying the chain rule gives us the Hessian: \[\mathcal{H}= d^2_q\mathcal{J}= \partial^2_q\mathcal{J}+ \partial^2_qG_{\boldsymbol{\lambda}} + (\partial_\mathbf{u}\partial_q\mathcal{J})d_q\mathbf{u}+ (\partial_\mathbf{u}\partial_qG_{\boldsymbol{\lambda}})d_q\mathbf{u}+ (\partial_{\boldsymbol{\lambda}}\partial_qG_{\boldsymbol{\lambda}})d_q\boldsymbol{\lambda}.\]

It is instructive to consider the action of this Hessian on two perturbations \(\delta q_1\) and \(\delta q_2\):

\[\label{Hessian95action} \begin{align} \mathcal{H}[\delta q_1][\delta q_2] = &(\partial^2_q\mathcal{J}+ \partial^2_qG_{\boldsymbol{\lambda}})[\delta q_1][\delta q_2] ~+ \\ &\langle~d_q\boldsymbol{\lambda}[\delta q_1],\partial_q\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\delta q_2]~\rangle ~+\\ &\langle~d_q\mathbf{u}[\delta q_1],(\partial_q\partial_\mathbf{u}\mathcal{J}+ \partial_q\partial_\mathbf{u}G_{\boldsymbol{\lambda}})[\delta q_2]~\rangle \end{align}\tag{29}\]

where we have used the fact that \((\partial_q\partial_\mathbf{u}\mathcal{F})^\ast = \partial_\mathbf{u}\partial_q\mathcal{F}\) for some functional \(\mathcal{F}\) (as well as other similar identities).
As with the first-order adjoint system, derivatives of the constraints at our disposal are employed to reduce this to an equation that doesn’t depend on Jacobians that would be prohibitive to compute. The new adjoint system will again constitute a system of equations whose solutions correspond to the Lagrange multipliers reqiured to do this.
To 29 , we add two indetically zero terms: the inner product of two Lagrange multipliers \(\boldsymbol{\mu}\) and \(\boldsymbol{\beta}\) with the gradients of the model and first-order adjoint equations applied to a perturbation \(\delta q_1\):

\[\begin{align} \mathcal{H}[\delta q_1][\delta q_2] = &(\partial^2_q\mathcal{J}+ \partial^2_qG_{\boldsymbol{\lambda}})[\delta q_1][\delta q_2] ~+ \\ &\langle~d_q\boldsymbol{\lambda}[\delta q_1],\partial_q\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\delta q_2]~\rangle ~+\\ &\langle~d_q\mathbf{u}[\delta q_1],(\partial_q\partial_\mathbf{u}\mathcal{J}+ \partial_q\partial_\mathbf{u}G_{\boldsymbol{\lambda}})[\delta q_2]~\rangle ~+\\ &\langle~\boldsymbol{\beta}, d_q G[\delta q_1]~\rangle ~+\\ &\langle~\boldsymbol{\mu}, d_q(\partial_{\mathbf{u}}\mathcal{J}+ \partial_{\mathbf{u}}G_{\boldsymbol{\lambda}})[\delta q_1]~\rangle. \end{align}\]

After expanding and rearranging, we find: \[\begin{align} \mathcal{H}[\delta q_1][\delta q_2] = &(\partial^2_q\mathcal{J}+ \partial^2_qG_{\boldsymbol{\lambda}})[\delta q_1][\delta q_2] ~+ \\ &\langle~\boldsymbol{\beta}, \partial_{\mathbf{u}}G[\delta q_1]~\rangle~+\\ &\langle~\boldsymbol{\mu}, (\partial_{q}\partial_{\mathbf{u}}\mathcal{J}+ \partial_{q}\partial_{\mathbf{u}}G_{\boldsymbol{\lambda}})[\delta q_1]~\rangle~+\\ &\langle~d_q\boldsymbol{\lambda}[\delta q_1],\partial_q\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\delta q_2] + \partial_{\mathbf{u}}\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\boldsymbol{\mu}]~\rangle ~+\\ &\langle~d_q\mathbf{u}[\delta q_1],(\partial_q\partial_\mathbf{u}\mathcal{J}+ \partial_q\partial_\mathbf{u}G_{\boldsymbol{\lambda}})[\delta q_2] + (\partial_{\mathbf{u}}G)^\ast[\boldsymbol{\beta}] + (\partial_{\mathbf{u}}^2G_{\boldsymbol{\lambda}}+ \partial_{\mathbf{u}}^2\mathcal{J})[\boldsymbol{\mu}]~\rangle. \end{align}\]

This shows that we can formulate a Hessian-vector product of the form:

\[\begin{align} \label{eq:svp95expansion} \mathcal{H}[\delta q] = (\partial^2_q\mathcal{J}+ \partial^2_qG_{\boldsymbol{\lambda}})[\delta q]+(\partial_\mathbf{u}\partial_q\mathcal{J}+ \partial_\mathbf{u}\partial_qG_{\boldsymbol{\lambda}})[\boldsymbol{\mu}] + (\partial_q G)^\ast[\boldsymbol{\beta}] \end{align}\tag{30}\]

if we choose the \(\boldsymbol{\mu}\) and \(\boldsymbol{\beta}\) that solve the equations:

\[\begin{align} &\partial_q\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\delta q_2] + \partial_{\mathbf{u}}\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\boldsymbol{\mu}] = 0\label{lniegkap}\\ &(\partial_q\partial_\mathbf{u}\mathcal{J}+ \partial_q\partial_\mathbf{u}G_{\boldsymbol{\lambda}})[\delta q_2] + (\partial_{\mathbf{u}}G)^\ast[\boldsymbol{\beta}] + (\partial_{\mathbf{u}}^2G_{\boldsymbol{\lambda}}+ \partial_{\mathbf{u}}^2\mathcal{J})[\boldsymbol{\mu}] = 0. \end{align}\tag{31}\]

These are the second-order adjoint equations.

9.2 Application to the SSA↩︎

To apply this to the shallow-stream approximation to the Cauchy momentum equations for ice sheet flow, there are 12 terms in the above equations that we need to evaluate. Choosing a cost functional that is the sum of separate functionals of \(\mathbf{u}\) and \(q\) means that terms involving \(\partial_\mathbf{u}\partial_q\mathcal{J}\) can be set to zero. As with the calculation of first-order derivatives, we assume a viscosity that is independent of \(\mathbf{u}\), so that \(\bar{\mu}\) doesn’t depend on \(\mathbf{u}\) or \(q\). This means that \(G\) is linear in \(\mathbf{u}\) and we can set \(\partial_\mathbf{u}^2 G_{\boldsymbol{\lambda}} =0\) as well. This leaves five expressions in the SOA system (now the “SOSA” system, after the afforementioned approximation) that need to be evaluated, and four in the equation for the Hessian-vector product. These are:

\[\begin{align} &\partial_q\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\delta q]\tag{32}\\ &\partial_{\mathbf{u}}\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\boldsymbol{\mu}]\tag{33}\\ &\partial_q\partial_{\mathbf{u}}G_{\boldsymbol{\lambda}}[\delta q]\tag{34}\\ &\partial_{\mathbf{u}}^2\mathcal{J}[\boldsymbol{\mu}]\tag{35}\\ &(\partial_{\mathbf{u}}G)^\ast[\boldsymbol{\beta}]\tag{36}\\ &\partial_q^2\mathcal{J}[\delta q]\tag{37}\\ &\partial_q^2G_{\boldsymbol{\lambda}}[\delta q]\tag{38}\\ &\partial_q\partial_\mathbf{u}G_{\boldsymbol{\lambda}}[\boldsymbol{\mu}]\tag{39}\\ &(\partial_q G)^\ast[\boldsymbol{\beta}]\tag{40}. \end{align}\]

We will deal with these in order, starting with the first 5 which define the SOSA system:

• Eq. 32 . \[\begin{align} \partial_{q}\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\delta q] &= \partial_{q}G[\delta q]\\ &= \lim_{\varepsilon\rightarrow0}\left\{\frac{1}{\varepsilon}\left(G(q+\varepsilon\delta q, \mathbf{u}) - G(q, \mathbf{u})\right)\right\}\\ &= \nabla\cdot (\varphi~\delta q~ H(\mathbf{u})) \end{align}\] where we have used that everything is smooth and linear.

• Eq. 33 :

  \[\begin{align} \partial_{\mathbf{u}}\partial_{\boldsymbol{\lambda}}G_{\boldsymbol{\lambda}}[\boldsymbol{\mu}] &= \partial_{\mathbf{u}}G[\boldsymbol{\mu}] = G(\boldsymbol{\mu})\\ &= \nabla\cdot(\varphi~ H(\boldsymbol{\mu}))-C\boldsymbol{\mu}. \end{align}\] This makes the crucial SOSA assumption that \(G\) is linear.

• Eq. 34 :

  Consider the action of \(\partial_{\mathbf{u}}G_{\boldsymbol{\lambda}}\) on \(\delta\mathbf{u}\). As \(G_{\boldsymbol{\lambda}}\) is linear, this is just \(G_{\boldsymbol{\lambda}}(\delta\mathbf{u})=\int\boldsymbol{\lambda}\cdot\left[\nabla\cdot(\varphi~\delta q~ H(\delta\mathbf{u})) - C\delta\mathbf{u}\right]~d\Omega\). Consider: \[\begin{align} \langle~ \partial_{q}G_{\boldsymbol{\lambda}}(\delta\mathbf{u}), \delta q ~\rangle &= \int\boldsymbol{\lambda}\cdot\nabla\cdot(\varphi~\delta q~ H(\delta\mathbf{u})) d\Omega\\ &= \int\delta\mathbf{u}\cdot\nabla\cdot(\varphi~\delta q~ H(\boldsymbol{\lambda})) d\Omega \label{final95but} \end{align}\tag{41}\] where the second line is once more due to the self-adjointness of the viscous part of \(G\). Eq. 41 shows the action of the mixed derivative \(\partial_{q}\partial_{\mathbf{u}}G_{\boldsymbol{\lambda}}\) on \((\delta q, \delta\mathbf{u})\). Hence the map \(\partial_{q}\partial_{\mathbf{u}}G_{\boldsymbol{\lambda}}[\delta q]\) can be seen to take the form: \[\partial_{q}\partial_{\mathbf{u}}G_{\boldsymbol{\lambda}}[\delta q] = \nabla\cdot(\varphi\delta q H(\boldsymbol{\lambda})).\]

• Eq 35 : We assume that, by construction, \(\partial_{\mathbf{u}}^2\mathcal{J}[\boldsymbol{\mu}]\) is not too difficult to write down.

• Eq. 36 : \[\begin{align} (\partial_{\mathbf{u}}G)^\ast[\boldsymbol{\beta}] &= \nabla\cdot(\varphi H(\boldsymbol{\beta})) - C\boldsymbol{\beta}. \end{align}\] This is the self-adjointness again.

With these in the bag, the second-order adjoint equations for the shallow-stream approximation to the Stokes equations are:

\[\begin{align} \label{soa95system} &\grad\cdot(\varphi\delta qH(\mathbf{u})) + \nabla\cdot(\varphi H(\boldsymbol{\mu})) - C\boldsymbol{\mu}= 0\\ &\grad\cdot(\varphi\delta qH(\boldsymbol{\lambda})) + \nabla\cdot(\varphi H(\boldsymbol{\beta})) - C\boldsymbol{\beta}+ \partial_{\mathbf{u}}^2\mathcal{J}[\boldsymbol{\mu}]= 0. \end{align}\tag{42}\]

For the case in which \(\mathcal{J}= \int(\mathbf{u}-\mathbf{u}_{\rm obs})^2~d\Omega\), the term \(\partial_{\mathbf{u}}^2\mathcal{J}[\boldsymbol{\mu}]\) is equal to \(2\boldsymbol{\mu}\).
Once again, the derivation of these equations has required the we set a litany of boundary terms to zero. Hence we require all variables obey the same boundary conditions as \(\mathbf{u}\)
Once these equations are solved, the Hessian-vector product can be evaluated according to eq. 30 . Let us derive 37 40 .

• Eq 37 . \(\partial_{q}^2\mathcal{J}[\delta q]\) will depend on the form of the regularisation used in the cost function.

• Eq 38 : \[\begin{align} \partial_{q}^2G_{\boldsymbol{\lambda}}[\delta q] &= \partial_{q}\langle\partial_{q}G_{\boldsymbol{\lambda}},\delta q\rangle\\ &=\partial_{q}\int\boldsymbol{\lambda}\cdot\nabla\cdot(\varphi~\delta q~ H(\mathbf{u}))~d\Omega\\ &= -\partial_q\int(\nabla\boldsymbol{\lambda}:(\varphi~\delta q~ H(\mathbf{u})))\\ &= -\nabla\boldsymbol{\lambda}:(\varphi~\delta q~ H(\mathbf{u})) \end{align}\]

• The same kind of reasoning gives: Eq. 39 : \[\begin{align} \partial_{\mathbf{u}}\partial_{q}G_{\boldsymbol{\lambda}}[\boldsymbol{\mu}] = -\nabla\boldsymbol{\lambda}:\varphi H(\boldsymbol{\mu}) \end{align}\]

• Eq. 40 : \[\begin{align} (\partial_{q}G)^\ast[\boldsymbol{\beta}] = -\nabla\boldsymbol{\beta}:\varphi H(\mathbf{u}). \end{align}\]

Hence, the Hessian-vector product can be calculated according to:

\[\begin{align} \label{eq:hvp95final} \mathcal{H}[\delta q] = \partial^2_q\mathcal{J}[\delta q] - \nabla\boldsymbol{\lambda}:(\varphi~\delta q~ H(\mathbf{u})) -\nabla\boldsymbol{\lambda}:\varphi H(\boldsymbol{\mu}) -\nabla\boldsymbol{\beta}:\varphi H(\mathbf{u}). \end{align}\tag{43}\]

Taking a linear rheology for the calculation of first and second order sensitivities has led to a system of equations for the adjoint variables \(\boldsymbol{\lambda}\), \(\boldsymbol{\mu}\), \(\boldsymbol{\beta}\) that all essentially have the same form, with different right-hand sides. This means little development is required to perform the calculations necessary for approximating Hessian-vector products once a forward SSA code has been set up. This makes the SOSA model derived here an appealing drop-in to many current ice sheet models, requiring only the use of existing stencils.

References↩︎