Preface

Numerical models of glaciers and ice sheets include parameters such as viscosity (often expressed using other parameters), density, thermal conductivity, heat capacity, geothermal heat flux and basal traction coefficients. (Even basal topography can be considered a model ‘parameter’.) Many of these parameters are difficult to measure directly. We use the phrase ‘data assimilation’ to refer to a method where more easily measured data, such as surface velocities, are used to estimate the values of such parameters. The idea behind data assimilation is to find the values of the parameter of interest (e.g. a frictional coefficient along the bed) that result in the best match between values calculated by the model (surface velocities, perhaps) and measured data. Thus, finding appropriate values for model parameters can be accomplished by solving an optimization problem, namely that of minimizing a measure of the difference between calculated and measured values.

In this paper ‘adjoint methods’ are described in the context of optimization problems. Application to an idealized data-assimilation problem from glaciology is described in the final section. Although applications using real glaciological data can be found in various published research papers (e.g. Reference Brinkerhoff, Meierbachtol, Johnson and HarperBrinkerhoff and others, 2011) no attempt is made to review them here.

A Typical Situation

Often the problem to be solved (whose solution is u = F(p)) is expressed as a system of n equations, f(u, p) = 0. For this problem we still have a function g (u), for which we want to find optimal values of the parameters p and thus need ∂g/∂ p = [∂g/∂p ₁ ∂g/∂p ₂ ∂g/∂p_m ]^T. In the preceding section the chain rule was used to show that ∂g/∂ p = ∂ u/∂p ^T ∂g/∂ u. To simplify the appearance of some of the equations that follow, we begin by taking the transpose of both sides of this equation and use the matrix identity (AB) ^T = B^T A^T to obtain

As before, the elements of ∂g/∂ u are obtained by differentiating g (u), but now the elements of ∂ u /∂ p must be determined from f(u, p) = 0. If we take the total derivative of both sides of the equation f_i (u, p) = 0 with respect to p_k we obtain

Doing this for every combination of the n elements of u and the m elements of p results in nm equations that can be represented compactly as ∂ f/ ∂ u ∂ u /∂ p + ∂ f/∂p = 0:

This equation can be solved for the ∂u_i/∂p_k factors

needed to obtain our goal

Here ∂g /∂u ^T is a 1 × n (row) matrix, ∂f/∂u is n × n and ∂f /∂p is n × m. To avoid multiplying the two large matrices we perform the multiplication in the order

With the product λ^T = ∂g /∂u^T ∂f/∂u⁻¹ found by solving the adjoint problem

Note that ∂ f /∂ u is a Jacobian matrix; once again obtaining the desired partial derivatives (∂g/∂ p) uses the transpose (i.e. adjoint) of a Jacobian matrix. Also note that this Jacobian matrix is the one used in Newton’s method to (approximately) solve f(u, p) = 0 (with fixed p) iteratively.

Example 2:

Suppose

The system of equations f ₁₍ u ₁, u ₂, p ₁, p ₂) = 0 and f ₂(u ₁, u ₂, p ₁, p ₂) = 0 has a solution, (u ₁, u ₂₎, that depends on the parameters p ₁ and p ₂. Changing the values of the parameters p ₁ and p ₂ causes the values of the variables u ₁ and u ₂ to change and thus the value of g (u ₁, u ₂) also changes. We want to find ∂g/∂ p ^T = [∂g/∂p ₁ ∂g/∂p ₂].

To calculate ∂g/∂ p ^T using the adjoint method we begin by finding the Jacobian matrix

along with

Next we write the adjoint problem (∂f /∂u)^T λ = ∂g /∂u

and solve it to find

The desired derivatives are now found from

To appreciate the meaning of these derivatives we consider a couple of specific value pairs for the parameters p ₁ and p ₂.

First, if p ₁ = 0 and p ₂ = 0, the equation f(u, p) = 0 has solution u = 0, so g( u) = 0. Since g = u ² ₁ + u ² ₂ ≥ 0 (always), g( u) = 0 must be a global minimum so we expect ∂g/∂ p = 0. Indeed, substituting u ₁ = 0 and u ₂ = 0 into our expression, confirm this,

To gain further appreciation of ∂g/∂ p, we consider a graphical representation of our example problem. If we plot solutions of f ₁(u ₁, u ₂, p ₁, p ₂₎ = u ₁ + u ₂ + p ₁ = 0 in the (u ₁, u ₂₎ plane we obtain a line with slope –1 and vertical intercept –p ₁. If we plot solutions of f ₂₍ u ₁, u ₂, p ₁, p ₂) = u ³ ₁ – u ₂ + p ₂ = 0 in the (u ₁, u ₂) plane we obtain a cubic curve that crosses the vertical axis at p ₂ with slope 0. The solution to the system of equations f(u, p) = 0 is represented by the point where the line and the cubic curve cross. The function g gives the square of the distance from the origin, so circles centered on the origin represent contours of constant g. Curves for p ₁ = –2, p ₂ = 0 and g = 2 are shown in Figure 1.

Fig. 1. Curves for the problem presented in example 2, with p ₁ = –2, p ₂ = 0 and g( u) = 2.

As can be seen from Figure 1, for this value of p the solution of f(u, p) = 0 is u = (1, 1). If we increase p ₁, the line with slope 1 will move down in the figure so the solution point will follow the cubic curve, decreasing its distance from the origin, so g decreases. If, however, we increase p ₂, the cubic curve will move up so the solution point will follow the line with slope –1 which is tangent to the (dashed) circle representing a contour of constant g. Thus we expect ∂g/∂p ₁ < 0 and ∂g=∂p ₂ = 0. Substituting the solution u ₁ = 1 and u ₂ = 1 into our expression gives ∂g/∂p ^T = [–2 0], confirming this. The reader is encouraged to use the graphical representation to find other parameter pairs where one of the components of ∂g/∂p is zero and use the expression for ∂g/∂p ^T to confirm their findings.

A Numerical Example

In this section we use finite differences to approximate the boundary-value problem

where p (x) = p ₀ + p ₁ x + p ₂ x ² (a second-order polynomial). After computing an approximate solution for particular values of the parameters c ₂, c ₁, c ₀, p ₀, p ₁, p ₂, a ₀ and a ₁, we use the adjoint method to calculate the rate of change of a function of the solution, g (u), as we vary these parameters. Two different functions are considered: (1) g is the value of u at the midpoint of the domain, and (2) g is the average value of u,

To apply the finite-difference method to our differential equation we divide the domain [0, 1] into n subintervals, each with length The endpoints of these subintervals are x ₀, x ₁, …, x_n , where x_k = k∆x. We denote approximate values of u at each of these points as u_k . The first derivative, u ^ʹ, at x_k is approximated by . The second derivative, u ^ʹʹ, at x_k is approximated by u_k ’ = Substituting these approximations into the differential equation gives

which can be rearranged to give

a linear algebraic equation that can be written for each of the interior points k = 1, 2, …, n –1. At the endpoints of the domain we assert u ₀ = a ₀ and u_n = a ₁. This collection of equations can be written compactly as Au = b, where A is a tridiagonal matrix, and u and b are column matrices. The matrix equation can be solved (e.g. using Gaussian elimination and back substitution) for u, giving an approximation to the solution of the boundary-value problem.

Following the procedure in the preceding paragraph for the (somewhat arbitrary) parameter values c ₂ = 1, c ₁ = –2, c ₀ = 1, p ₀ = 1, p ₁ = 1, p ₂ = –5, a ₀ = 0 and a ₁ = 0 generates the approximation shown in Figure 2.

Fig. 2. Solution to the boundary-value problem u″ – 2u′ + u = 1 + x – 5x ² for 0 < x < 1, and u( 0) = 0, u( 1) = 0. The solid curve is the exact solution; the black circles are an approximate solution found using finite differences with n = 20.

We can use the adjoint method to find the rate of change of a function of u with respect to the parameters c ₂, c ₁, c ₀, p ₀, p ₁, p ₂, a ₀ and a ₁. We consider first the simple function where an even value for n has been chosen.

To utilize the procedure described in the preceding section, the system of equations represented by the matrix equation Au = b is written in the form f(u, p) = 0, by subtracting the right-hand sides from both sides of each equation. Here p = [c ₂ c ₁ c ₀ p ₀ p ₁ p ₂ a ₀ a ₁]. To find ∂g/∂ p we need the three matrices ∂ f/ @ u, ∂ f/ @ p and ∂g/∂ u. Since ∂f/ ∂ u and ∂ f/ ∂ p require two columns to be displayed clearly, they are shown in Figure 3. For g (u) = u_{n/ 2},

Fig. 3. Matrices used to calculate ∂g/∂ p for the numerical example involving a boundary-value problem.

After solving the adjoint problem ∂ f/ ∂ u ^T λ = ∂g/∂ u for λ we can calculate ∂g/∂ p ^T = –λ^T∂ f/∂p. The results for the example problem whose solution is shown in Figure 2 are given in Table 1.

The values of the partial derivatives found using the adjoint method were verified by perturbing the parameters one at a time. For each parameter (c ₂, c ₁, c ₀, p ₀, p ₁, p ₂, a ₀ and a ₁) an additional finite-difference approximation was calculated with that parameter increased by 0.01, giving new values of u. Calling these new values ũ we can approximate the partial derivative of g with respect to the perturbed parameter by . As seen in Table 1, these approximations are very close to the values found using the adjoint method.

To calculate the partial derivatives of we use (the composite) Simpson’s rule (which requires an even value for n) to approximate this integral,

from which we have

Again, ∂g/∂p ^T = –λ^T ∂f/∂p, where λ is the solution of the adjoint problem ∂f/∂u ^Tλ = ∂g/∂u. Results are given in Table 1.

A Data-Assimilation Example from Glaciology

Mathematical models of glaciers typically include parameters which are difficult to estimate directly from field measurements. For example, some type of frictional coefficient is required at the base of the glacier where the ice rests upon the earth. The basal traction (described mathematically using this frictional coefficient) affects the velocity of the ice throughout the glacier. Thus, it is plausible to use the more easily measured velocity at the upper surface of the glacier to determine the frictional coefficient. In this concluding section we present an example problem, using a simplified geometry to demonstrate the use of the adjoint method to accomplish this.

The example problem is related to a benchmark problem, ‘Experiment D’ of the Ice Sheet Model Intercomparison Project for Higher-Order Models (ISMIP-HOM) (Reference PattynPattyn and others, 2008). In ISMIP-HOM Experiment D, a sheet of ice with uniform thickness of 1000 m rests on a plane inclined at an angle α = 0: 1°. A domain of length 20 km with periodic boundary conditions at each end is considered. A numerical model is used to calculate the velocity of the ice, which varies not only throughout the thickness of the ice but also along the length of the domain, due to variations in basal traction. There is no variation in the horizontal direction transverse to the inclination.

In ISMIP-HOM Experiment D, the basal traction is prescribed and ice velocities (and pressure) are calculated. We call this a ‘forward problem’. Here we are interested in a related ‘inverse problem’: given the surface velocity, find the basal traction.

The model used is a finite-element approximationFootnote ^* of Stokes’ equationsFootnote ^†

along with conservation of mass,

In these equations, x and y are the coordinate directions, x along the incline and y perpendicular to it (almost vertical). The components of velocity in the x – and y –directions are given by u and v, respectively. Pressure is represented by the variable p, while ρ = 910 kg m^–3 is the density of ice and g = 9: 81 m s⁻² is the acceleration due to gravity. The viscosity, μ, varies according to Glen’s flow lawFootnote ^‡

with the effective strain rate

And flow parameters A = 10⁻¹⁶ Pa^-na⁻¹ and n = 3.

Stress-free boundary conditions apply at the upper surface

While at the base v = 0 and basal traction is modeled using the boundary condition

In ISMIP-HOM Experiment D

where L is the length of the domain (20 km for the problem presented here). This frictional coefficient is plotted as a heavy black curve in the lower left panel of Figure 4. The upper left panel of Figure 4 shows the velocity component, u, at the (top) surface of the ice for a solution of the ‘forward’ problem (using β ² = 1000 + 1000 sin 2πx/L)) plotted as a heavy black line.

Fig. 4. Solution of an inverse problem. The component of the ice velocity in the direction parallel to the bed (u; m a⁻¹) at the (top) surface of the ice is shown in the upper left panel. This velocity evolves towards the desired solution (the heavy black curve) as the basal friction (β ², plotted in the lower left panel) is changed. The upper right panel shows the error, g ðu Þ, decreasing as the coefficients, p_k , in the trigonometric expansion of β ² (shown in the lower right panel) change.

Now, for the inverse problem. Suppose we are given the surface velocity component, u, plotted as a heavy black curve in the upper left panel of Figure 4. We hope to use this information to discover the basal friction coefficient, β ² (plotted as a heavy black curve in the lower left panel of Fig. 4). To do so we pose the problem as an optimization problem: find β ² such that the surface velocity component, u, found by solving the ‘forward’ problem minimizes the error given by

Where u_d is the desired surface velocity. To discretize the function β² a trigonometric expansionFootnote ^§ is used,

A variety of algorithms for solving optimization problems exist (Reference PressPress, 2007). Here we choose the Broyden–Fletcher– Goldfarb–Shano (BFGS) algorithm (Reference PressPress, 2007, p. 521–524), which uses partial derivatives of the objective function, g(u), with respect to the parameters, p_k , that can be varied. These partial derivatives (∂g/∂p_k ) are exactly what the adjoint method described in the previous sections can provide. Starting with a ‘guess’ for β ², an approximate solution to the inverse problem is obtained by repeating the following steps:

1. 1. Solve the ‘forward’ problem with the current estimate of β².

2. 2. Evaluate the error, g(u).

3. 3. Evaluate the partial derivatives, ∂g/∂p_k , using the adjoint method.

4. 4. Calculate an improved estimate of β² (i.e. the p_k ) according to the BFGS algorithm.

The results of applying this process are shown in Figure 4. The initial guess was β² = 1000; i.e. p ₀ = 1000 and p_k = 0 for 1 ≤ k ≤ N where, here, N = 30. The gray curves in the left panels of Figure 4 show approximations of u and β ² as the process progressed. The original guess of² ¼ 1000 (represented by a horizontal line in the lower left panel of Fig. 4) resulted in the bottommost horizontal line in the upper left panel of Figure 4. As the estimate of² improves, the curves approach the desired values plotted as heavy black curves. The upper right panel of Figure 4 shows the error, g (u), decreasing throughout the process. Further improvement can be made by allowing more iterations of the algorithm. The lower right panel of Figure 4 shows the evolution of the coefficient (p_k ) values. The upper curve represents p ₀ which starts at 1000, dips down for a while but returns to a value very near 1000. The curve that increases from zero to near 1000 represents p ₁. The remainder of the curves, which remain near zero, represent p ₂ through p ₃₀. In the optimal solution these values are all zero, but the largest, p ₉, was 36.9 when the process was terminated. Again, further improvement can be made by allowing more iterations.