Parameter sensitivity analysis of dynamic ice sheet models-Numerical computations

The friction coefficient and the base topography of a stationary and a dynamic ice sheet are perturbed in two models for the ice: the full Stokes equations and the shallow shelf approximation. The sensitivity to the perturbations of the velocity and the height at the surface is quantified by solving the adjoint equations of the stress and the height equations providing weights for the perturbed data. The adjoint equations are solved numerically and the sensitivity is computed in several examples in two dimensions. Comparisons are made with analytical solutions to simplified problems.


Introduction
The result of isothermal simulations of large ice sheets depends on the ice model, the topography, and the parametrization of the conditions at the base of the ice. The models are systems of partial differential equations (PDEs) for the velocity, pressure, and height of the ice. The topography and the friction model with its parameters determine the horizontal velocity and the height at the ice surface in the computations. In the inverse problem, the parameters at the base are inferred from data at the surface by solving adjoint equations and minimizing the difference between given data and simulated results. In this paper, we estimate the sensitivity of the surface observations to changes in the basal conditions by solving the adjoint equations to the full Stokes (FS) equations and the shallow shelf (or shelfy stream) approximation (SSA), see Greve and Blatter (2009);MacAyeal (1989). The advantage of solving the adjoint equations in a variational control method is that the effect of many perturbations of the parameters at the bottom is obtained for one observation at one point of the surface at a certain time point. If there are many observations and only one perturbation, then it is more efficient to compute the sensitivity by solving the forward model PDEs twice in a direct method, firstly with the unperturbed parameters, secondly with the perturbed parameters, and then take the difference between the solutions. The direct method has the advantage that there is no need to implement a solver for the adjoint equations.
Most methods for inversion of ice surface data to compute parameters in the models at the ice base rely on a solution of the adjoint stress equation with a given fixed geometry of the ice as in MacAyeal (1993); Petra et al. (2012). The time dependent height equation for the moving upper surface is not included in the inversion. The stationary basal friction coefficients have been derived from satellite data in this way for many glaciers and continental ice sheets using velocity data in e.g. Gillet-Chaulet et al. (2016); Isaac et al. (2015); Schannwell et al. (2019); Sergienko and Hindmarsh (2013). The sensitivity to changes at the base increases closer to the grounding line in the coastal regions in Durand et al. (2011). The base topography is inferred from height data in van Pelt et al. (2013) without solving the adjoint equations. The conditions between the ice and the bedrock vary in time and sometimes the friction parameter varies several orders of magnitude in a decade in Jay- Allemand et al. (2011).
In addition, there are variations on seasonal and diurnal time scales with examples in Schoof (2010); Shannon et al. (2013); Vallot et al. (2017). Other time dependent forces are considered in Seddik et al. (2019). The effect of a seasonal variation of the lubrication at the base of the ice is studied in Shannon et al. (2013) for the Greenland ice sheet by solving the FS and other high order equations. Fast temporal variations in the meltwater under the ice drive the ice flow in the analysis in Schoof (2010). The spatial and temporal variations of the basal conditions are inferred from satellite data in Larour et al. (2014) with an inverse method for SSA and automatic differentiation. Based on observations, the conclusion in Sole et al. (2011) is also that the annual change of the water drainage under the ice affects the sliding and the acceleration and deceleration of the ice. Here, we solve the adjoint equations to both the stress equation and the time dependent height equation in FS and SSA to examine how the dynamics of the models change the sensitivity to the base parameters. The adjoint equations are derived and analytical solutions are found to simplified equations in a companion paper by Cheng and Lötstedt (2019).
The forward advection equation for the height and the stress equations for the velocity for FS are here solved numerically in two dimensions (2D) with Elmer/Ice ; Gillet-Chaulet et al. (2012)). The solver of the adjoint stress equation in Elmer/Ice is amended by the adjoint height equation. The forward and adjoint SSA equations are solved in 2D by a finite difference method. The perturbations are observed in the velocity and the height at certain points in space and time. Comparisons are made for steady state and time dependent problems between a direct calculation of the change at the ice surface and using the control technique with the adjoint solution. Simplified adjoint stress equations have been proposed and used in Martin and Monnier (2014); Morlighem et al. (2013); Mosbeux et al. (2016). The sensitivity in the SSA model is evaluated here for such simplifications in the adjoint SSA equations. The numerical solutions are also compared to the analytical formulas in Cheng and Lötstedt (2019). There is a transfer matrix between the perturbations in the parameters at the base and the observations at the surface. The properties of this matrix are evaluated to see which combinations of perturbations and observations that are well and ill-conditioned. In an ill-conditioned problem, the sensitivity is low at the surface to perturbations at the base. This matrix can be used to quantify the uncertainty in the ice flow due to uncertainties in the model parameters, see e.g. Bulthuis et al. (2019);Schlegel et al. (2018);Smith (2014).
The ice equations and the corresponding adjoint equations for FS and SSA are given in Sect. 2. The computed sensitivities are compared for the direct method and the control method in Sect. 3 for steady state and time dependent problems in 2D. The ice configuration is taken from the MISMIP benchmark project in Pattyn et al. (2012). The results are discussed and conclusions are drawn in Sections 4 and 5. Formulas from Cheng and Lötstedt (2019) are found in Appendix A.
Vectors and matrices are written in bold as a and A. The operations ⊗, :, and on vectors a and c, matrices A and C, and four index tensors A are defined by The norm of a vector a is defined by a = (a · a) 1/2 .
The equations of two ice models and their adjoint equations are stated in this section. The FS equations are considered to be an accurate model of ice sheets and the SSA equations are an approximation of the FS equations suitable e.g. for fast flowing ice on the ground and ice floating on water, see Greve and Blatter (2009).

Full Stokes equations
The FS equations are a system of PDEs for the velocity of the ice u(x, t) = (u 1 , u 2 , u 3 ) T , the pressure p(x, t), and the height h(x, y, t) with the coordinates x = (x, y, z) and time t. There is a stress equation satisfied by u and p and an advection equation for h. The adjoint equation of the stress equation is derived in Petra et al. (2012) and the adjoint equations of the stress and the height equations are found in Cheng and Lötstedt (2019). The sensitivity of observations of the velocity and the height of the ice surface is derived for perturbations in the friction coefficient at the ice base.
The domain of the ice is Ω with boundary Γ in three dimensions (3D). The boundary consists of the ice surface at the upper boundary Γ s , the lower boundary at the ice base Γ b and Γ w , and the vertical, lateral boundaries Γ u and Γ d where Γ u is the upstream boundary with n · u ≤ 0 and Γ d is the downstream boundary with n · u > 0. The normal of Γ pointing outward is denoted by n. The projection of Γ s and Γ b on the horizontal x − y plane is ω and the projections of Γ u and Γ d are γ u and γ d , respectively. The z coordinate of the grounded base Γ b is the topography and the bathymetry b(x, y). The grounding line γ GL separates Γ b on ω from Γ w floating on water with a moving z-coordinate z b (x, y, t). Formal definitions of these domains are Let I be the identity matrix. The projection of a vector on the tangential plane of Γ b is denoted by T = I − n ⊗ n as in Petra et al. (2012). In 2D, x = (x, z) T , ω = [0, L], γ u = 0, and γ d = L.

Forward equations
The definitions of the strain rate D and the viscosity η of the ice are The trace of D 2 is trD 2 and the rate factor A depends on the temperature of the ice, here assumed to be constant in isothermal flow. The material constant n > 0 is given in Glen's flow law. Then the stress tensor is Let ρ be the density of the ice, g be the gravitational acceleration and a be the accumulation/ablation rate on the surface Γ s .
The notation is simplified with the slope vectors h = (h x , h y , −1) T in 3D and h = (h x , −1) T in 2D. A subscript x, y, z, or t on a variable denotes a partial derivative such that e.g. h x = ∂h/∂x. Then the forward FS equations for h, u, and p are The initial data for h are h 0 (x) and h γ (x, t) is specified on the inflow boundary γ u . The expression Cf (Tu) defines the friction law with variable coefficient C(x, t) and a function f (·) of the projected velocity Tu, e.g. as in Weertman (1957) where The Dirichlet boundary conditions of u on Γ u and Γ d are set to be u u and u d .

Adjoint equations
We observe a quantity at the surface Γ s when t ∈ [0, T ]. For example, if the ice is in the steady state and F (u) = u 1 δ(x − x * ) with the Dirac delta δ then the observation is the x component of u at The adjoint equations depend on the first variations F u and F h of F (u, h) with respect to u and h. In the first example above, , 0, 0) T and F h = 0 and in the second example F u = 0 and F h = δ(x − x * ).
The adjoint FS equations form a system of PDEs for the adjoint height ψ, the adjoint velocity v, and the adjoint pressure q.
There is an advection equation for ψ and an adjoint stress equation for v and q such that where the adjoint viscosity, adjoint stress, and linearized friction law in Eq. (8) are according to Petra et al. (2012) The tensor I with four indices ijkl is 1 when i = j = k = l and 0 otherwise.
The perturbation of the observation in Eq. (7) with respect to a perturbation in the friction coefficient C is involving the tangential projections of the forward and adjoint velocities Tu and Tv at the grounded ice base Γ b . This expression is derived in Cheng and Lötstedt (2019) and Petra et al. (2012) via the perturbation of the Lagrangian of the system of equations and evaluating it at the forward and adjoint solutions.
Only perturbations in C are considered here for the FS model. Via the Lagrangian, the result of perturbations δb in the topography can be derived but the complexity of the adjoint Eq. (8) would increase considerably.

Shallow shelf approximation
In the shallow shelf approximation of the FS equations, the velocity is constant in the vertical direction and the pressure is given by the cryostatic approximation (Greve and Blatter (2009);MacAyeal (1989)). The sensitivity of observations of the velocity at the surface and the height to perturbations in friction coefficients and the base topography is quantified for the SSA model.

Forward equations
It is sufficient to solve for the horizontal velocity u = (u 1 , u 2 ) T when x = (x, y) ∈ ω thus simplifying the 3D FS problem Eq. (5) considerably. The viscosity in the SSA is where B(u) = D(u) + ∇ · u I. The stress tensor ς(u) in SSA is defined by Let n be the outward normal vector of the boundary γ, t the tangential vector such that n · t = 0, and H = h − b the thickness of the ice. The friction law is defined as in the FS case in Eq. (6) where the basal velocity is replaced by the horizontal velocity since the vertical variation is neglected in SSA. Under the floating ice shelf on Γ w , C = 0 in the friction law.
The ice dynamics system is where u in ≤ 0 and u out > 0 are the inflow and outflow normal velocities on γ u and γ d of the boundary γ = γ u ∪γ d . The friction on the lateral side of the ice γ = γ g ∪ γ w depends on the tangential velocity t · u there. The friction law C γ f γ (t · u) on γ g is not necessarily the same as Cf (u) on ω.
The structure of the SSA system Eq. (13) is similar to the FS equations in Eq. (5). However, the velocity u is not divergence free in SSA and B = D due to the cryostatic approximation.

Adjoint equations
The adjoint SSA equations are derived in Cheng and Lötstedt (2019) cf.η andσ in Eq. (9). The adjoint SSA equations are With a Weertman friction law Eq. (6), the terms F ω and F γ in the adjoint basal friction and the lateral friction in Eq. (15) are The friction coefficients on the base and the lateral sides are perturbed by δC and δC γ and the topography is perturbed by δb in the SSA model. Then the perturbation δF in the observation F in Eq. (7) is (Cheng and Lötstedt (2019))

Forward and adjoint SSA in 2D
In the 2D model, u 2 = 0, derivatives with respect to y vanish, and the lateral friction force is neglected, C γ = 0.
Assume that u > 0 and u x > 0. There is an inflow of ice with speed u L to the left and a calving rate u c at x = L. The viscosity in Eq. (11) is simplified to η = 2A −1/n u ν x . The friction term is Cf (u)u = Cu m with the Weertman law in Eq. (6). The adjoint variables v and ψ satisfy the adjoint equations in 2D Perturbations δb and δC in the topography and the friction coefficient propagate to the surface as in Eq. (16)

Discretized relations in 2D
In order to simplify the notation, only a 2D steady state problem for the SSA model is considered here but the analysis is applicable to 3D steady state problems as well as time-dependent problems with the FS or SSA models.
The time independent perturbation of F in Eq. (19) for the steady state solution is rewritten with F u = δ(x − x * ) and weights w ub and w uC The weights w ub and w uC in Eq. (20) depend on both x * and x. When h is observed the perturbation is where the weights w hb and w hC have the same form as in Eq. (20) but with different ψ and v.
The relation is discretized by observing u at equidistant (20) is computed by the trapezoidal rule to have or in matrix form with the matrix elements In the same manner, there are matrices W hb and W hC connecting δh with δb and δC The sensitivity of u to changes in b and C on ω is given by the singular value decomposition (SVD) of W ub and W uC (Golub and Loan (1989)) defined by where U ub and U uC are of size M ×M and V ub and V uC are of size N ×N . They are orthogonal matrices, e.g. U T ub U ub = I. The diagonal matrices Σ ub and Σ uC are of size M × N with non-negative singular values σ ubi and σ uCi in the diagonals ordered from large to small for increasing i = 1, 2, ..., min(M, N ).
Consider a case with δb = 0, the perturbation is simplified to δu = W uC δC. If M = N and the smallest singular value σ uCN = min i σ uCi is positive then If M > N with more observations of δu i than discrete δC j , then δC for a given δu can be computed in the least squares sense by minimizing δu − W uC δC also with the solution where Σ −1 uC is the generalized inverse of Σ uC of dimension N × M with elements σ −1 uCi on the diagonal and 0 elsewhere. The relation between δu and δC is well behaved in Eq. (25) and Eq. (26) if all the singular values σ uCi are of similar size, but if some of them are much smaller than the other ones with σ Ci σ C1 , i = J, J + 1, . . . , min(M, N ), then the relation is ill-conditioned. A large perturbation in C may then result in a hardly visible perturbation at the surface and a small observed perturbation in u may correspond to a large perturbation at the base. The same conclusions apply to W ub and σ ubi in the relation between δu and δb and to the sensitivity matrices W hb and W hC when The transfer functions in Gudmundsson (2003) between perturbations in b and C at the base and the observations u and h at the top are determined by linearization and Fourier transformation in a slab geometry. The transfer function for different wave numbers corresponds to the singular values in our analysis.

Results
In the numerical experiments we use a 2D constant downward-sloping bed with an ice profile from the MISMIP benchmark project in Pattyn et al. (2012). The bedrock elevation in meters is given as The initial configuration of the ice is a steady state solution achieved by the FS model using Elmer/Ice  with A = 1.38 × 10 −24 s −1 Pa −3 with a grounding line position at x GL = 1.053 × 10 6 m shown in Fig. 1. The Weertman type friction law in Eq. (6) in the forward problem has the exponent m = 1/3 and a constant friction coefficient C 0 = 7.624 × 10 6 m −1/3 s 1/3 Pa. The remaining physical parameters are given in Table 1 Without losing the generality in the friction law and to investigate the relation between the basal velocity and the stress, the friction law exponent in the adjoint problem is assumed to be m = 1 and the coefficient is calculated from the forward The multiplier ψ only acts as the amplitude of the external force on Γ s and h is an approximate normal vector pointing inward on Γ s in the adjoint FS equation Eq. (8). The size of ψh is several orders of magnitude smaller than 1, the coefficient in front of Consequently, in the u 1 -response case, the adjoint solution v is mainly influenced by the observation function F u . However, in the h-response case with F u = 0, the adjoint solution v is determined by ψh and the solution would be v = 0 if we did not solve the adjoint advection equation for ψ. The amplitude of the perturbation at the surface depends on the wavelength λ of the perturbation at the base. The shorter λ is, the smaller the amplitude is. Introduce a stationary perturbation δC(x) = C 0 cos(2π(x − x * )/λ) with a constant C 0 and a small 1. Then the change in the steady state solution u 1 at the surface is according to Eq. (10) The same relation holds for δh(x * ) but with a different v. Let be a measure of the width of the weight function for the steady state in Fig. 2 which is about 10 5 . When λ is large compared to then which is a constant value for long λ, and the perturbation can be observed at the surface. If the wavelength of the basal perturbation is short compared to , then it is damped before it reaches the surface and the effect of δC on u 1 and h is small.

SSA
The same MISMIP benchmark experiment as in Sect The analytical approximation of u is poor to the right of x GL for the floating ice in Fig. 5 but we are only interested in the solution on the ground. The reason for the error in the analytical solution of u is that H is assumed to be constant for x > x GL .
The analytical solution for H catches the fast decrease when x approaches x GL from the left. Another solution for x > x GL is found in Greve and Blatter (2009) assuming that the thickness depends linearly on x. largest and smallest singular values of Σ uC are 10 −4 and 10 −12 with a large quotient σ uC1 /σ uCN and the span of the singular values of Σ hC is from 10 −4 to 10 −8 (which is better).
The singular values of the sensitivity matrices W ub and W hb in Fig. 9 are in the interval 10 −4 to 10 −7 from large to small.
They are better conditioned than the sensitivity matrices for C. In particular, Σ hb (in pink-red) in the h-response case has the lowest variation of the singular values. The inverse problem of solving for the topography b from the surface elevation h in the steady state setup is a well-posed problem compared to inferring C from u.
The same perturbation on C(x) as in Fig. 4 is imposed in the SSA simulations. The perturbed solutions after one year and 15,000 years (which is close to a steady state) are computed with the forward equations and then the reference solutions at the steady state without any perturbation are subtracted. This difference is compared to the perturbations obtained with the adjoint equations as in Fig. 4. In the one year perturbation experiment in Fig. 10, the transient weight functions in the upper panels in and 7. The weights can be approximated by −θ(x, t)δ (x − x * ) for some θ > 0. Then the surface response will be where δC jumps discontinuously at x = 0.9 × 10 6 and x = 1.0 × 10 6 . The same phenomenon is found for FS in Fig. 4 with an explanation in Fig. 2.
The perturbations δu and δh in the steady state in Fig. 12 have discontinuous derivatives δu x and δh x where δC has jumps. This is explained by the integral terms in (A3) and (A5). The discontinuities in the upper panel of Fig. 13 are caused by the jumps in δb at 0.9 × 10 6 and 1.0 × 10 6 and the first term in (A3). The jumps in δh in the lower panel of Fig. 13 are due to the first term in (A5).
All the predicted solutions from the adjoint SSA are in good agreement with the forward perturbation. The singular values of the transfer matrices corresponding to the two simplifications are shown in Fig. 15 where the two transfer matrices are denoted by W uC for the system coupling ψ and v and by W uC for the adjoint equation without ψ with a fixed H. The singular values in Σ uC are similar to those in Σ uC in Fig. 9 since the influence of the adjoint viscosity on the system is almost negligible. The transfer matrix W uC has a better conditioning than W uC , although it is still worse than the best cases in Fig. 9. This implies that the inversion of steady state SSA without the height coupling may be an ill-posed problem.
Regularization is necessary penalising oscillatory behavior at the base as in Gagliardini et al. (2013); Petra et al. (2012).

Discussion
A few issues are discussed here related to the control method for estimating the parameter sensitivity. We solve the FS adjoint problem only one step backward in time to verify the numerical method due to limitations of the current framework of Elmer/Ice. It is possible but more complicated and expensive to solve the adjoint problem numerically for a large number of time steps K. This requires storing all the forward solutions (u i , p i , h i ), i = 1, 2, . . . , K, to be able to compute the adjoint solutions (v i , q i , ψ i ), i = K, K − 1, . . . , 1, which may be prohibitive in 3D. Since the data to be stored in the SSA model is one dimension lower, we are able to solve the adjoint problem backward in time for any number of K.
However, for a fair comparison, we show the results for one time step with SSA in this paper. The transfer relation W uC between small perturbations of the friction coefficient C at the ice base and the perturbation of the horizontal velocity u at the ice surface is given by Eq. (23) with δb = 0. The singular values of W uC in Fig. 9 tell how sensitive u is to changes in C. The transfer relation also describes how the uncertainty in C is propagated to uncertainty in the velocity at the surface and how uncertainty δu in measurements of u appear as uncertainty δC in C Eq. (26), see Smith (2014).
The transfer relation is computed by solving the forward problem once and then the adjoint problem for each one of the M observations. An alternative would be to solve the forward equations first for the unperturbed solution and then perturb C by δC j and solve the forward equations again N times and subtract to find the relation between δu and δC j . It is usually more expensive to solve the nonlinear forward equations than the linear adjoint equations. Suppose that the computational work to solve the forward problem is W F and the adjoint problem is W A . If the forward and adjoint equations are in similar form, such as the FS or SSA problem, and solving the nonlinear forward problem requires k iterations where every nonlinear iteration has the same computational cost as solving the linear adjoint problem, then W A /W F ≈ 1/k. The quotient between the work to determine the transfer relation involving the adjoint equations and the work only based on the forward equation is solve N + 1 forward problems to compute W uC . In the inverse problem to find C given observations of u, h, the functions F u and F h are smooth and M = 1 in the iterative procedure to compute C. Solving the adjoint equations is then always favorable.

Conclusions
The part of u is observed and has little influence on δu. On the contrary, if h is observed then there is an important effect of ψ on δh in FS and SSA. The magnitudes of ψ are different depending on whether u or h is observed. Simplifications of the SSA adjoint in the steady state by using the forward viscosity or ignoring the adjoint height equation have minor consequences for the predictions of u with a perturbed C in Fig. 14. The sensitivity to perturbations δb and δC is quantified for steady state and time dependent problems with the FS and SSA models. It increases as the observation point x * approaches the grounding line. This is explained by analytical expressions for SSA where the sensitivity is inversely proportional to the ice thickness H(x * ). The closer we are to the grounding line the higher the requirements are on the resolution of the topography and the friction coefficient to obtain accurate solutions of u and h there.
A weight is local if its extension in space is close to the observation point. The weights on δC at the ice base are local for the steady state and time dependent FS model. They are also local for the time dependent SSA model and the transfer from δb to δu and δh in the steady state. The sensitivity of δu and δh in the steady state of SSA depends on δC from a larger domain. It is difficult to observe a perturbation δC with a short wavelength on u and h. In the example in Fig. 3, a spatial perturbation wavelength λ = 2 × 10 4 m (about 10H) in C is damped by 0.2 in δh and 0.02 in δu compared to a wavelength λ > 10 5 where there is no damping due to λ. The transfer matrices from δb and δC to δu and δh are examined by the singular value decomposition. If the quotient between the largest and the smallest singular values of the matrix is large then it is ill-conditioned and if it is small (but ≥ 1) then the problem is well-conditioned. In an ill-conditioned problem, some perturbations at the base will be barely visible at the surface and a small perturbation at the top may correspond to a large perturbation at the bottom. In a well-conditioned problem, all perturbations at the base have a measurable effect at the surface. The ranking of the conditioning of the transfers in Fig. 9 from the best to the worst is 1. δb → δh, 2. δb → δu, 3. δC → δh, 4. δC → δu.
In the past, the coupling between δu and δC is most frequently used for inference of C from velocity data but height data could improve the robustness of the inference.