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. A transfer matrix couples the perturbations at the base with the perturbations at the top. Comparisons are made with analytical solutions to simplified problems. The sensitivity to perturbations depends on their wavelengths and the distance to the grounding line. A perturbation in the topography has a direct effect at the ice surface above it, while a change in the friction coefficient is less visible there.

1 Introduction

The output 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 boundary conditions of the PDEs are given by the topography and the friction model with its parameters. Of particular interest in the simulation of ice is the horizontal velocity and the height at the ice surface. 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 and 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.

We are interested in the effect of perturbations of the topography and the slipperiness at the ice base on the velocity of the ice at the surface and its height. By solving the adjoint equations, we quantify the sensitivity to perturbations close to the grounding line and of different wavelengths. The sensitivity at the upper surface to perturbations in the basal topography and friction is different, and the separation of the two contributions appears to be difficult. The transfer functions between the perturbations at the base and the surface observations are more or less well behaved. A related problem is to infer the basal geometry and friction coefficients from observational data by inversion using the adjoint solution.

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) for SSA and in Petra et al. (2012) for FS, where 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), and Sergienko and Hindmarsh (2013).

The conditions between the ice and the bedrock vary in time, and sometimes the friction parameter varies several orders of magnitude in a decade (see, e.g., 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), and 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 in the water drainage under the ice affects the sliding and the acceleration and deceleration of the ice. Transient data are included in Goldberg et al. (2015) to find time-dependent basal parameters by inversion, where the sensitivity is determined by automatic differentiation. The results differ if the time evolution of the equations is taken into account or not. The shallow ice approximation (SIA) is the ice model in Monnier and des Boscs (2017) to determine the basal properties with time-dependent surface data. 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 influence of the dynamics of the basal conditions is different on the velocity and the height observations.

The forward advection equation for the height and the stress equations for the velocity for FS are here solved numerically in two dimensions (2-D) with Elmer/Ice (Gagliardini et al., 2013; 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 for a vertical ice in 2-D 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), and Mosbeux et al. (2016). The sensitivity in the SSA model is evaluated in this paper for such simplifications in the adjoint SSA equations. The sensitivity in the numerical solutions is also compared to the analytical formulas in Cheng and Lötstedt (2019). It is observed in Durand et al. (2011) that the sensitivity to changes at the base increases closer to the grounding line in the coastal regions. The basal topography is inferred from the height data in van Pelt et al. (2013) without solving the adjoint equations. The reason for the increased sensitivity and why the height method works are explained by our analytical solutions to the adjoint SSA equations.

There is a transfer matrix between the perturbations in the parameters at the base and the observations at the surface. Analytical expressions for time-dependent transfer functions for FS and SSA are derived in Gudmundsson (2003, 2008) by linearizing, freezing coefficients, and applying Fourier analysis and the Laplace transform. The properties of the transfer matrix are evaluated here to see which combinations of perturbations and observations are well- and ill-conditioned. In an ill-conditioned problem, the sensitivity at the surface to perturbations at the base is low. This matrix can also 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). Perturbations at the ice base with short wavelength are propagated to the surface with a weaker effect on the height and velocity compared to long wavelengths in Gudmundsson (2003, 2008). These are the conclusions in calculations with FS in Kyrke-Smith et al. (2018), where it is difficult to separate the contribution from the friction and the bed topography from each other. These effects are confirmed in our analysis.

The structure of the paper is as follows. The ice equations and the corresponding adjoint equations for FS and SSA are presented in Sect. 2. The computed sensitivities are compared between the direct method and the control method in Sect. 3 for steady-state and time-dependent problems in 2-D. The ice configuration is taken from the MISMIP benchmark project in Pattyn et al. (2012). The results are discussed and conclusions are drawn in Sects. 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 $\otimes ,:,$ and ⋆ on vectors a and c, matrices A and C, and four index tensors 𝒜 are defined by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }(\mathit{a}\otimes \mathit{c}{)}_{ij}=a_i{c}_j,\mathit{a}:\mathit{c}=\mathit{a}\cdot \mathit{c}=\sum _i{a}_i{c}_i,\\ {\displaystyle }& {\displaystyle }(\mathbf{A}\otimes \mathbf{C}{)}_{ijkl}=A_{ij}{C}_{kl},\mathbf{A}:\mathbf{C}=\sum _{ij}{A}_{ij}{C}_{ij},\\ \text{(1)}& {\displaystyle }& {\displaystyle }(\mathcal{A}\star \mathrm{C}{)}_{ij}=\sum _{kl}{\mathcal{A}}_{ijkl}{C}_{kl}.\end{array}$$

The definition of a norm of a vector a is $\Vert \mathit{a}\Vert =(\mathit{a}\cdot \mathit{a}{)}^{\mathrm{1}/\mathrm{2}}$.

2 Ice models

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).

2.1 Full Stokes equations

The FS equations are a system of PDEs for the velocity of the ice $\mathit{u}(\mathit{x},t)=(u_{\mathrm{1}},u_{\mathrm{2}},u_{\mathrm{3}}{)}^T$, the pressure p(x,t), and the height $h(x,y,t)$ with the coordinates $\mathit{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 (3-D). 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 and lateral boundaries Γ_u and Γ_d, where Γ_u is the upstream boundary with $\mathit{n}\cdot \mathit{u}\le \mathrm{0}$ and Γ_d is the downstream boundary with $\mathit{n}\cdot \mathit{u}>\mathrm{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

$$\begin{array}{ll}{\displaystyle }\mathrm{\Omega }& {\displaystyle }=\mathit{\left\{}\mathit{x}\right|(x,y)\in \mathit{\omega },b(x,y)\le z\le h(x,y,t)\mathit{\}},\\ {\displaystyle }{\mathrm{\Gamma }}_s& {\displaystyle }=\mathit{\left\{}\mathit{x}\right|(x,y)\in \mathit{\omega },z=h(x,y,t)\mathit{\}},\\ {\displaystyle }{\mathrm{\Gamma }}_b& {\displaystyle }=\mathit{\left\{}\mathit{x}\right|(x,y)\in \mathit{\omega },z=b(x,y),x<x_{\text{GL}}\left(y\right)\mathit{\}},\\ {\displaystyle }{\mathrm{\Gamma }}_w& {\displaystyle }=\mathit{\left\{}\mathit{x}\right|(x,y)\in \mathit{\omega },z=z_b(x,y,t),x>x_{\text{GL}}\left(y\right)\mathit{\}},\\ {\displaystyle }{\mathrm{\Gamma }}_u& {\displaystyle }=\mathit{\left\{}\mathit{x}\right|(x,y)\in {\mathit{\gamma }}_u,b(x,y)\le z\le h(x,y,t)\mathit{\}},\\ \text{(2)}& {\displaystyle }{\mathrm{\Gamma }}_d& {\displaystyle }=\mathit{\left\{}\mathit{x}\right|(x,y)\in {\mathit{\gamma }}_d,b(x,y)\le z\le h(x,y,t)\mathit{\}}.\end{array}$$

Let I be the identity matrix. The projection of a vector on the tangential plane of Γ_b is denoted by $\mathbf{T}=\mathbf{I}-\mathit{n}\otimes \mathit{n}$ as in Petra et al. (2012). In 2-D, $\mathit{x}=(x,z{)}^T,\mathit{\omega }=[\mathrm{0},L],{\mathit{\gamma }}_u=\mathrm{0},$ and γ_d=L.

2.1.1 Forward equations

The definitions of the strain rate D and the viscosity η of the ice are

$$\begin{array}{}\text{(3)}& \begin{array}{r}\mathbf{D}=\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^T),\mathit{\eta }\left(\mathit{u}\right)=\frac{\mathrm{1}}{\mathrm{2}}{A}^{-\frac{\mathrm{1}}{n}}\left(\text{tr}{\mathbf{D}}^{\mathrm{2}}\right(\mathit{u}){)}^{\mathit{\nu }},\mathit{\nu }=\frac{\mathrm{1}-n}{\mathrm{2}n}.\end{array}\end{array}$$

The trace of D² is trD² 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

$$\begin{array}{}\text{(4)}& \mathit{\sigma }(\mathit{u},p)=\mathrm{2}\mathit{\eta }\mathbf{D}\left(\mathit{u}\right)-p\mathbf{I}.\end{array}$$

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 $\mathit{h}=(h_x,h_y,-\mathrm{1}{)}^T$ in 3-D and $\mathit{h}=(h_x,-\mathrm{1}{)}^T$ in 2-D. A subscript $x,y,z,$ or t on a variable denotes a partial derivative such that $h_x=\partial h/\partial x$. Then the forward FS equations for $h,\mathit{u},$ and p are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{h}_t+\mathit{h}\cdot \mathit{u}=a,\text{on}{\mathrm{\Gamma }}_s,\\ {\displaystyle }& {\displaystyle }h(\mathit{x},\mathrm{0})=h_{\mathrm{0}}\left(\mathit{x}\right),\mathit{x}\in \mathit{\omega },h(\mathit{x},t)=h_{\mathit{\gamma }}(\mathit{x},t),\mathit{x}\in {\mathit{\gamma }}_u,\\ {\displaystyle }& {\displaystyle }-\mathrm{\nabla }\cdot \mathit{\sigma }(\mathit{u},p)=-\mathrm{\nabla }\cdot \left(\mathrm{2}\mathit{\eta }\right(\mathit{u}\left)\mathbf{D}\right(\mathit{u}\left)\right)+\mathrm{\nabla }p=\mathit{\rho }\mathit{g},\\ {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0},\text{in}\mathrm{\Omega }\left(t\right),\mathit{\sigma }\mathit{n}=\mathbf{0},\text{on}{\mathrm{\Gamma }}_s,\\ \text{(5)}& {\displaystyle }& {\displaystyle }\mathbf{T}\mathit{\sigma }\mathit{n}=-Cf\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u},\mathit{n}\cdot \mathit{u}=\mathrm{0},\text{on}{\mathrm{\Gamma }}_b.\end{array}$$

The initial data for h are h₀(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

$$\begin{array}{}\text{(6)}& f\left(\mathit{u}\right)=\Vert \mathit{u}{\Vert }^{m-\mathrm{1}},m>\mathrm{0}.\end{array}$$

The Dirichlet boundary conditions of u on Γ_u and Γ_d are set to be u_u and u_d.

2.1.2 Adjoint equations

We observe a quantity

$$\begin{array}{}\text{(7)}& \mathcal{F}=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_s}{\int }F(\mathit{u},h)\mathrm{d}\mathit{x}\mathrm{d}t\end{array}$$

at the surface Γ_s when $t\in [\mathrm{0},T]$. For example, if the ice is in the steady state and $F\left(\mathit{u}\right)=u_{\mathrm{1}}\mathit{\delta }(\mathit{x}-{\mathit{x}}_{\ast })$ with the Dirac delta (δ), then the observation is the x component of u at x_∗:

$${\displaystyle }{\displaystyle }\mathcal{F}=\underset{{\mathrm{\Gamma }}_s}{\int }F\left(\mathit{u}\right)\mathrm{d}\mathit{x}=u_{\mathrm{1}}\left({\mathit{x}}_{\ast }\right).$$

If $F\left(h\right)=h\mathit{\delta }(\mathit{x}-{\mathit{x}}_{\ast })$ then the height is observed as

$${\displaystyle }{\displaystyle }\mathcal{F}=\underset{{\mathrm{\Gamma }}_s}{\int }F\left(h\right)\mathrm{d}\mathit{x}=h\left({\mathit{x}}_{\ast }\right).$$

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, $F_u=\left(\mathit{\delta }\right(\mathit{x}-{\mathit{x}}_{\ast }),\mathrm{0},\mathrm{0}{)}^T$ and F_h=0, and in the second example F_u=0 and $F_h=\mathit{\delta }(\mathit{x}-{\mathit{x}}_{\ast })$.

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_t+\mathrm{\nabla }\cdot \left(\mathit{u}\mathit{\psi }\right)-\mathit{h}\cdot {\mathit{u}}_z\mathit{\psi }=F_h+F_u\cdot {\mathit{u}}_z,\text{on}{\mathrm{\Gamma }}_s,\\ {\displaystyle }& {\displaystyle }\mathit{\psi }(\mathit{x},T)=\mathrm{0},\mathit{\psi }(\mathit{x},t)=\mathrm{0},\text{on}{\mathrm{\Gamma }}_d,\\ {\displaystyle }& {\displaystyle }-\mathrm{\nabla }\cdot \stackrel{\mathrm{̃}}{\mathit{\sigma }}(\mathit{v},q)=-\mathrm{\nabla }\cdot \left(\mathrm{2}\stackrel{\mathrm{̃}}{\mathit{\eta }}\right(\mathit{u})\star \mathbf{D}(\mathit{v}\left)\right)+\mathrm{\nabla }q=\mathbf{0},\\ {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathit{v}=\mathrm{0},\text{in}\mathrm{\Omega }\left(t\right),\\ {\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\sigma }}(\mathit{v},q)\mathit{n}=-(F_{\mathit{u}}+\mathit{\psi }\mathit{h}),\text{on}{\mathrm{\Gamma }}_s,\\ {\displaystyle }& {\displaystyle }\mathbf{T}\stackrel{\mathrm{̃}}{\mathit{\sigma }}(\mathit{v},q)\mathit{n}=-Cf\left(\mathbf{T}\mathit{u}\right)\left(\mathbf{I}+{\mathbf{F}}_b\left(\mathbf{T}\mathit{u}\right)\right)\mathbf{T}\mathit{v},\text{on}{\mathrm{\Gamma }}_b,\\ \text{(8)}& {\displaystyle }& {\displaystyle }\mathit{n}\cdot \mathit{v}=\mathrm{0},\text{on}{\mathrm{\Gamma }}_b,\end{array}$$

where the adjoint viscosity, adjoint stress, and linearized friction law in Eq. (8) are according to Petra et al. (2012):

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{u}\right)& {\displaystyle }=\mathit{\eta }\left(\mathit{u}\right)\left(\mathcal{I}+{\displaystyle \frac{\mathrm{1}-n}{n\mathbf{D}\left(\mathit{u}\right):\mathbf{D}\left(\mathit{u}\right)}}\mathbf{D}\left(\mathit{u}\right)\otimes \mathbf{D}\left(\mathit{u}\right)\right),\\ {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\sigma }}(\mathit{v},q)& {\displaystyle }=\mathrm{2}\stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{u}\right)\star \mathbf{D}\left(\mathit{v}\right)-q\mathbf{I},\\ \text{(9)}& {\displaystyle }{\mathbf{F}}_b\left(\mathbf{T}\mathit{u}\right)& {\displaystyle }={\displaystyle \frac{m-\mathrm{1}}{\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{u}}}\left(\mathbf{T}\mathit{u}\right)\otimes \left(\mathbf{T}\mathit{u}\right).\end{array}$$

The tensor ℐ 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

$$\begin{array}{}\text{(10)}& \mathit{\delta }\mathcal{F}=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathrm{\Gamma }}_b}{\int }f\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{v}\mathit{\delta }C\mathrm{d}\mathit{x}\mathrm{d}t,\end{array}$$

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 by 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.

2.2 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.

2.2.1 Forward equations

It is sufficient to solve for the horizontal velocity $\mathit{u}=(u_{\mathrm{1}},u_{\mathrm{2}}{)}^T$ when $\mathit{x}=(x,y)\in \mathit{\omega }$, thus simplifying the 3-D FS problem Eq. (5) considerably. The viscosity in the SSA is

$$\begin{array}{ll}{\displaystyle }\mathit{\eta }\left(\mathit{u}\right)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{-\frac{\mathrm{1}}{n}}{\left(u_{\mathrm{1}x}^{\mathrm{2}}+u_{\mathrm{2}y}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}(u_{\mathrm{1}y}+u_{\mathrm{2}x}{)}^{\mathrm{2}}+u_{\mathrm{1}x}{u}_{\mathrm{2}y}\right)}^{\mathit{\nu }}\\ \text{(11)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{-\frac{\mathrm{1}}{n}}{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathbf{B}:\mathbf{D}\right)}^{\mathit{\nu }},\end{array}$$

where $\mathbf{B}\left(\mathit{u}\right)=\mathbf{D}\left(\mathit{u}\right)+\mathrm{\nabla }\cdot \mathit{u}\mathbf{I}$. The stress tensor ς(u) in SSA is defined by

$$\begin{array}{}\text{(12)}& \mathit{\varsigma }\left(\mathit{u}\right)=\mathrm{2}H\mathit{\eta }\mathbf{B}\left(\mathit{u}\right).\end{array}$$

Let n be the outward normal vector of the boundary γ, t the tangential vector such that $\mathit{n}\cdot \mathit{t}=\mathrm{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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{h}_t+\mathrm{\nabla }\cdot \left(\mathit{u}H\right)=a,\mathrm{0}\le t\le T,\mathit{x}\in \mathit{\omega },\\ {\displaystyle }& {\displaystyle }h(\mathit{x},\mathrm{0})=h_{\mathrm{0}}\left(\mathit{x}\right),\mathit{x}\in \mathit{\omega },h(\mathit{x},t)=h_{\mathit{\gamma }}(\mathit{x},t),\mathit{x}\in {\mathit{\gamma }}_u,\\ {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathit{\varsigma }-Cf\left(\mathit{u}\right)\mathit{u}=\mathit{\rho }gH\mathrm{\nabla }h,\mathit{x}\in \mathit{\omega },\\ {\displaystyle }& {\displaystyle }\mathit{n}\cdot \mathit{u}(\mathit{x},t)=u_{\text{in}}(\mathit{x},t),\mathit{x}\in {\mathit{\gamma }}_u,\\ {\displaystyle }& {\displaystyle }\mathit{n}\cdot \mathit{u}(\mathit{x},t)=u_{\text{out}}(\mathit{x},t),\mathit{x}\in {\mathit{\gamma }}_d,\\ {\displaystyle }& {\displaystyle }\mathit{t}\cdot \mathit{\varsigma }\mathit{n}=-C_{\mathit{\gamma }}{f}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})\mathit{t}\cdot \mathit{u},\mathit{x}\in {\mathit{\gamma }}_g,\\ \text{(13)}& {\displaystyle }& {\displaystyle }\mathit{t}\cdot \mathit{\varsigma }\mathit{n}=\mathrm{0},\mathit{x}\in {\mathit{\gamma }}_w,\end{array}$$

where u_in≤0 and u_out>0 are the inflow and outflow normal velocities on γ_u and γ_d of the boundary $\mathit{\gamma }={\mathit{\gamma }}_u\cup {\mathit{\gamma }}_d$. The friction on the lateral side of the ice $\mathit{\gamma }={\mathit{\gamma }}_g\cup {\mathit{\gamma }}_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.

2.2.2 Adjoint equations

The adjoint SSA equations are derived in Cheng and Lötstedt (2019) as in Sect. 2.1.2 by forming the Lagrangian and partial integration using the forward equations and the boundary conditions in Eq. (13). The adjoint viscosity $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ and adjoint stress $\stackrel{\mathrm{̃}}{\mathit{\varsigma }}$ are defined by

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{u}\right)& {\displaystyle }=\mathit{\eta }\left(\mathit{u}\right)\left(\mathcal{I}+{\displaystyle \frac{\mathrm{1}-n}{n\mathbf{B}\left(\mathit{u}\right):\mathbf{D}\left(\mathit{u}\right)}}\mathbf{B}\left(\mathit{u}\right)\otimes \mathbf{D}\left(\mathit{u}\right)\right),\\ \text{(14)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\varsigma }}\left(\mathit{v}\right)& {\displaystyle }=\mathrm{2}H\stackrel{\mathrm{̃}}{\mathit{\eta }}\left(\mathit{u}\right)\star \mathbf{B}\left(\mathit{v}\right);\end{array}$$

cf. $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ and $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ in Eq. (9). The adjoint SSA equations are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_t+\mathit{u}\cdot \mathrm{\nabla }\mathit{\psi }+\mathrm{2}\mathit{\eta }\mathbf{B}\left(\mathit{u}\right):\mathbf{D}\left(\mathit{v}\right)-\mathit{\rho }gH\mathrm{\nabla }\cdot \mathit{v}\\ {\displaystyle }& {\displaystyle }+\mathit{\rho }g\mathit{v}\cdot \mathrm{\nabla }b=F_h,\text{in}\mathit{\omega },\\ {\displaystyle }& {\displaystyle }\mathit{\psi }(\mathit{x},T)=\mathrm{0},\text{in}\mathit{\omega },\mathit{\psi }(\mathit{x},t)=\mathrm{0},\text{on}{\mathit{\gamma }}_w,\\ {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \stackrel{\mathrm{̃}}{\mathit{\varsigma }}\left(\mathit{v}\right)-Cf\left(\mathit{u}\right)(\mathbf{I}+{\mathbf{F}}_{\mathit{\omega }}(\mathit{u}\left)\right)\mathit{v}-H\mathrm{\nabla }\mathit{\psi }=-F_u,\text{in}\mathit{\omega },\\ {\displaystyle }& {\displaystyle }\mathit{t}\cdot \stackrel{\mathrm{̃}}{\mathit{\varsigma }}\left(\mathit{v}\right)\mathit{n}=-C_{\mathit{\gamma }}{f}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})(\mathrm{1}+F_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u}\left)\right)\mathit{t}\cdot \mathit{v},\text{on}{\mathit{\gamma }}_g,\\ {\displaystyle }& {\displaystyle }\mathit{t}\cdot \stackrel{\mathrm{̃}}{\mathit{\varsigma }}\left(\mathit{v}\right)\mathit{n}=\mathrm{0},\text{on}{\mathit{\gamma }}_w,\\ \text{(15)}& {\displaystyle }& {\displaystyle }\mathit{n}\cdot \mathit{v}=\mathrm{0},\text{on}\mathit{\gamma }.\end{array}$$

Compared to Eq. (8), the advection equation depends on v, and the influence of ψ in the stress equation is different in Eq. (15). 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

$${\displaystyle }{\displaystyle }{\mathbf{F}}_{\mathit{\omega }}\left(\mathit{u}\right)={\displaystyle \frac{m-\mathrm{1}}{\mathit{u}\cdot \mathit{u}}}\mathit{u}\otimes \mathit{u},F_{\mathit{\gamma }}=m-\mathrm{1}.$$

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 δℱ in the observation ℱ in Eq. (7) is (Cheng and Lötstedt, 2019)

$$\begin{array}{ll}{\displaystyle }\mathit{\delta }\mathcal{F}=& {\displaystyle }\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathit{\omega }}{\int }\left(\mathrm{2}\mathit{\eta }\mathbf{B}\right(\mathit{u}):\mathbf{D}(\mathit{v})+\mathit{\rho }g\mathit{v}\cdot \mathrm{\nabla }h+\mathrm{\nabla }\mathit{\psi }\cdot \mathit{u})\mathit{\delta }b\\ {\displaystyle }& {\displaystyle }-f\left(\mathit{u}\right)\mathit{u}\cdot \mathit{v}\mathit{\delta }C\mathrm{d}\mathit{x}\mathrm{d}t\\ \text{(16)}& {\displaystyle }& {\displaystyle }-\underset{\mathrm{0}}{\overset{T}{\int }}\underset{{\mathit{\gamma }}_g}{\int }{f}_{\mathit{\gamma }}(\mathit{t}\cdot \mathit{u})\mathit{t}\cdot \mathit{u}\mathit{t}\cdot \mathit{v}\mathit{\delta }{C}_{\mathit{\gamma }}\mathrm{d}s\mathrm{d}t.\end{array}$$

2.2.3 Forward and adjoint SSA in 2-D

In the 2-D model, u₂=0, derivatives with respect to y vanish, and the lateral friction force is neglected, C_γ=0. The ice domains are the grounded and floating parts ${\mathrm{\Gamma }}_b=[\mathrm{0},x_{\text{GL}}]$ and ${\mathrm{\Gamma }}_w=(x_{\text{GL}},L]$, where x_GL is the position of the grounding line. The friction coefficient C is positive on Γ_b and C=0 on Γ_w. The forward and adjoint equations in 2-D are derived from Eqs. (13) and (15) by letting H and u₁ be independent of y and taking u₂=0. The notation is simplified if we let u=u₁ and v=v₁. The forward equations follow from Eq. (13):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{h}_t+(uH{)}_x=a,\mathrm{0}\le t\le T,\mathrm{0}\le x\le L,\\ {\displaystyle }& {\displaystyle }h(x,\mathrm{0})=h_{\mathrm{0}}\left(x\right),h(\mathrm{0},t)=h_L\left(t\right),\\ {\displaystyle }& {\displaystyle }(H\mathit{\eta }{u}_x{)}_x-Cf(u)u-\mathit{\rho }gH{h}_x=\mathrm{0},\mathrm{0}\le x\le L,\\ \text{(17)}& {\displaystyle }& {\displaystyle }u(\mathrm{0},t)=u_L\left(t\right),u(L,t)=u_c\left(t\right).\end{array}$$

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 $\mathit{\eta }=\mathrm{2}{A}^{-\mathrm{1}/n}{u}_x^{\mathit{\nu }}$. 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 2-D:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_t+u{\mathit{\psi }}_x+(\mathit{\eta }{u}_x-\mathit{\rho }gH)v_x+\mathit{\rho }g{b}_{x}v=F_h,\\ {\displaystyle }& {\displaystyle }\mathrm{0}\le t\le T,\mathrm{0}\le x\le L,\\ {\displaystyle }& {\displaystyle }{\left({\displaystyle \frac{\mathrm{1}}{n}}H\mathit{\eta }{v}_x\right)}_x-Cmf\left(u\right)v-H{\mathit{\psi }}_x=-F_u,\\ \text{(18)}& {\displaystyle }& {\displaystyle }\mathit{\psi }(x,T)=\mathrm{0},\mathit{\psi }(L,t)=\mathrm{0},v(\mathrm{0},t)=\mathrm{0},v(L,t)=\mathrm{0},\end{array}$$

obtained from Eqs. (14) and (15) or derived from Eq. (17) with equal result.

Perturbations δb and δC in the topography and the friction coefficient propagate to the surface as in Eq. (16):

$$\begin{array}{}\text{(19)}& \mathit{\delta }\mathcal{F}=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{0}}{\overset{L}{\int }}({\mathit{\psi }}_{x}u+v_x\mathit{\eta }{u}_x+v\mathit{\rho }g{h}_x)\mathit{\delta }b-vf\left(u\right)u\mathit{\delta }C\mathrm{d}x\mathrm{d}t.\end{array}$$

2.2.4 Discretized relations in 2-D

In order to simplify the notation, only a 2-D steady-state problem for the SSA model is considered here, but the analysis is applicable to 3-D steady-state problems as well as time-dependent problems with the FS or SSA models.

The time-independent perturbation of ℱ in Eq. (19) for the steady-state solution is rewritten with $F_u=\mathit{\delta }(x-x_{\ast })$ and weights w_ub and w_uC:

$$\begin{array}{ll}{\displaystyle }\mathit{\delta }u\left(x_{\ast }\right)& {\displaystyle }={\displaystyle }\mathit{\delta }\mathcal{F}=\underset{\mathrm{0}}{\overset{L}{\int }}{w}_{ub}\mathit{\delta }b+w_{uC}\mathit{\delta }C\mathrm{d}x,\\ {\displaystyle }{w}_{ub}(x_{\ast },x)& {\displaystyle }={\displaystyle }{\mathit{\psi }}_{x}u+v_x\mathit{\eta }{u}_x+v\mathit{\rho }g{h}_x,w_{uC}(x_{\ast },x)\\ \text{(20)}& {\displaystyle }& {\displaystyle }=-vf\left(u\right)u.\end{array}$$

The weights w_ub and w_uC in Eq. (20) depend on both x_∗ and x. When h is observed the perturbation is

$$\begin{array}{}\text{(21)}& \mathit{\delta }h\left(x_{\ast }\right)={\displaystyle }\underset{\mathrm{0}}{\overset{L}{\int }}{w}_{hb}\mathit{\delta }b+w_{hC}\mathit{\delta }C\mathrm{d}x,\end{array}$$

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 $x_{\ast i},i=\mathrm{1},\mathrm{2},\mathrm{\dots },M,$ with $x_{\ast ,i+\mathrm{1}}-x_{\ast i}=\mathrm{\Delta }{x}_{\ast }$, and perturbing b and C at $x_j,j=\mathrm{1},\mathrm{2},\mathrm{\dots },N,$ with $x_{j+\mathrm{1}}-x_j=\mathrm{\Delta }x$. The integral in Eq. (20) is computed by the trapezoidal rule to have

$$\begin{array}{ll}{\displaystyle }\mathit{\delta }u\left(x_{\ast i}\right)=& {\displaystyle }\sum _{j=\mathrm{1}}^N{\mathit{\mu }}_j\left(w_{ub}\right(x_{\ast i},x_j\left)\mathit{\delta }b\right(x_j)\\ {\displaystyle }& {\displaystyle }+w_{uC}(x_{\ast i},x_j)\mathit{\delta }C\left(x_j\right))\mathrm{\Delta }x,\\ \text{(22)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{1}}=& {\displaystyle }\mathrm{0.5},{\mathit{\mu }}_j=\mathrm{1},j=\mathrm{2},\mathrm{3},\mathrm{\dots },N-\mathrm{1},{\mathit{\mu }}_N=\mathrm{0.5},\end{array}$$

or in matrix form

$$\begin{array}{}\text{(23)}& \mathit{\delta }\mathit{u}={\mathbf{W}}_{ub}\mathit{\delta }\mathit{b}+{\mathbf{W}}_{uC}\mathit{\delta }\mathit{C},\end{array}$$

with the matrix elements

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{W}_{ubij}={\mathit{\mu }}_j{w}_{ub}(x_{\ast i},x_j),W_{uCij}={\mathit{\mu }}_j{w}_{uC}(x_{\ast i},x_j),\\ {\displaystyle }& {\displaystyle }i=\mathrm{1},\mathrm{2},\mathrm{\dots },M,j=\mathrm{1},\mathrm{2},\mathrm{\dots },N.\end{array}$$

In the same manner, there are matrices W_hb and W_hC connecting δh with δb and δC:

$$\begin{array}{}\text{(24)}& \mathit{\delta }h={\mathbf{W}}_{hb}\mathit{\delta }\mathit{b}+{\mathbf{W}}_{hC}\mathit{\delta }\mathit{C}.\end{array}$$

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 Van Loan, 1989) defined by

$${\displaystyle }{\displaystyle }{\mathbf{W}}_{ub}={\mathbf{U}}_{ub}{\mathbf{\Sigma }}_{ub}{\mathbf{V}}_{ub}^T,{\mathbf{W}}_{uC}={\mathbf{U}}_{uC}{\mathbf{\Sigma }}_{uC}{\mathbf{V}}_{uC}^T,$$

where U_ub and U_uC are of size M×M and V_ub and V_uC are of size N×N. A summary of the properties of the SVD is found in the Appendix.

Consider the case when δb=0 in Eq. (23). The relation between δu and δC is well behaved in Eq. (23) if all the singular values σ_uCi of W_uC are of similar size, but if some of them are much smaller than the other ones with ${\mathit{\sigma }}_{Ci}\ll {\mathit{\sigma }}_{C\mathrm{1}},i=J,J+\mathrm{1},\mathrm{\dots },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 $F_h=\mathit{\delta }(x-x_{\ast })$.

The transfer functions in Gudmundsson (2003, 2008) 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.

2.2.5 Relation to the inverse problem

The sensitivity problem and the inverse problem are related. Assume that there are M observations of the velocity u_obs at the surface of the ice at x_i and we want to derive the corresponding friction coefficient C at j locations. With C we observe u at the top at the same coordinates. Then we seek a correction δC of C at N points such that u+δu approaches u_obs. Using Eq. (23),

$$\begin{array}{}\text{(25)}& \mathit{u}-{\mathit{u}}_{\text{obs}}=\mathit{\delta }\mathit{u}={\mathbf{W}}_{uC}\mathit{\delta }\mathit{C},\end{array}$$

and δC is chosen such that $\Vert \mathit{u}-{\mathit{u}}_{\text{obs}}\Vert$ is minimized. This problem is a linear-least-squares problem. Expressed with the SVD and the generalized inverse ${\mathbf{\Sigma }}_{uC}^{-\mathrm{1}}$, the solution is

$$\begin{array}{}\text{(26)}& \mathit{\delta }\mathit{C}={\mathbf{V}}_{uC}{\mathbf{\Sigma }}_{uC}^{-\mathrm{1}}{\mathbf{U}}_{uC}^T(\mathit{u}-{\mathit{u}}_{\text{obs}});\end{array}$$

see Golub and Van Loan (1989) and the Appendix. The solution can be improved iteratively with updates of C and u by ${\mathit{C}}_{k+\mathrm{1}}={\mathit{C}}_k+\mathit{\delta }\mathit{C}$ and ${\mathit{u}}_{k+\mathrm{1}}={\mathit{u}}_k+\mathit{\delta }\mathit{u}$, computing a new W_uC and so on.

The relation between the transfer matrix and the inversion problem is illustrated by Eq. (26), but a more efficient optimization method is based on the gradient of the objective function $\Vert \mathit{u}-{\mathit{u}}_{\text{obs}}\Vert$. It is the standard method for inversion in, e.g., Gillet-Chaulet et al. (2016), Isaac et al. (2015), and Petra et al. (2012), and the gradient is computed using the adjoint solution with $F\left(\mathit{u}\right)=\Vert \mathit{u}-{\mathit{u}}_{\text{obs}}{\Vert }^{\mathrm{2}}$.

3 Results

In the numerical experiments we use a 2-D 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

$$\begin{array}{}\text{(27)}& b\left(x\right)=\mathrm{720}-\mathrm{778.5}\times {\displaystyle \frac{x}{\mathrm{7.5}\times {\mathrm{10}}^{\mathrm{5}}}}.\end{array}$$

Figure 1 The initial ice geometry with height h (blue), ice base b (orange), and ocean bathymetry (black). The domains in Eq. (2) are the ice domain Ω between the blue and orange curves, the upper surface Γ_s in blue, the lower boundary on the bedrock Γ_b and on water Γ_w in orange, Γ_u at x=0, and Γ_d at $x=L=\mathrm{1.6}\times {\mathrm{10}}^{\mathrm{6}}$ m.

The initial configuration of the ice is a steady-state solution achieved by the FS model using Elmer/Ice (Gagliardini et al., 2013) with $A=\mathrm{1.38}\times {\mathrm{10}}^{-\mathrm{24}}$ ${\mathrm{s}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{3}}$, with a grounding line position at $x_{\text{GL}}=\mathrm{1.053}\times {\mathrm{10}}^{\mathrm{6}}$ m shown in Fig. 1. The Weertman-type friction law in Eq. (6) in the forward problem has the exponent $m=\mathrm{1}/\mathrm{3}$ and a constant friction coefficient $C_{\mathrm{0}}=\mathrm{7.624}\times {\mathrm{10}}^{\mathrm{6}}$ ${\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}{\mathrm{s}}^{\mathrm{1}/\mathrm{3}}\mathrm{Pa}$. The remaining physical parameters are given in Table 1.

Table 1 The physical parameters of the ice.

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 steady-state solution by $C\left(\mathit{x}\right)=C_{\mathrm{0}}\Vert \mathit{u}{\Vert }^{-\mathrm{2}/\mathrm{3}}$. The resulting friction law becomes Cf(u)=C(x), which can be viewed as a linearization of the friction law at the steady state.

3.1 Full Stokes model

A vertically extruded mesh is constructed for the given geometry with mesh size Δx=1 km yielding equidistant nodes in the horizontal direction. The number of vertical layers is set to 20 in the whole domain. Only the grounded ice is considered in the adjoint problem, and Dirichlet boundary conditions on u are used for the lateral boundaries Γ_d and Γ_u at the grounding line x=x_GL and the ice divide x=0.

The forward and adjoint FS problems are solved using the finite element code Elmer/Ice (Gagliardini et al., 2013) with a P1–P1 quadrilateral element and Galerkin least-squares stabilization for the Stokes equation and a bubble stabilization (Baiocchi et al., 1993) for the adjoint advection equation. The feature to solve the adjoint time-dependent equations has been added to Elmer/Ice. The Dirac delta is approximated by a linear basis function with the amplitude 1∕Δx.

The time-stepping scheme for the forward and adjoint transient problems is the implicit Euler method with a constant time step Δt=1 year. The adjoint equation is solved backward in time from the final time t=T to t=0. The steady state of the adjoint equations is computed by neglecting the time derivative term in the adjoint surface equation Eq. (8) and solving the corresponding linear system of equations for ψ and v.

Both transient and steady-state simulations are run with pointwise observations of the horizontal velocity u₁ and surface elevation h at different x_∗ positions on the top surface. The time interval for the transient solutions is [0,1], covered by one forward time step Δt from 0 to 1 and one backward time step from 1 to 0.

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, which is equivalent to the coefficient in front of $\mathit{\delta }(x-x_{\ast })$ in F_u. Consequently, in the u₁-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 ψ.

Figure 2 Comparison of the weights Tu⋅Tv in Eq. (10) for perturbations δC at different observation points $x_{\ast }=\mathrm{0.25}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.5}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$, and 0.9×10⁶ (blue, orange, green, and pink). (a, b) Transient simulations; (c, d) steady states. (a, c) w_uC with pointwise u response; (b, d) w_hC with pointwise h response.

The adjoint solutions v₁ at Γ_b of all the four cases are concentrated at the observation points. The vertical component v₂ shares the same feature as v₁ due to the boundary condition $\mathit{n}\cdot \mathit{v}=\mathrm{0}$ on Γ_b. Therefore, the weights Tu⋅Tv in Fig. 2 are also confined to the neighborhood of x_∗. The negative weights obtained in the u₁-response cases imply that an increase in the basal friction coefficient results in a decrease in the surface velocity. The amplitude of the weights grows rapidly toward the grounding line in all four cases in the figure. In fact, the contribution of the weight function to the observed variables u₁ can be viewed as a convolution of the perturbation in C(x) with a narrow Gaussian $w_{uC}(x_{\ast },x)$ in Eq. (20) after a proper scaling in the left panels of Fig. 2.

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 $\mathit{\delta }C\left(x\right)=\mathit{ϵ}{C}_{\mathrm{0}}\cos \left(\mathrm{2}\mathit{\pi }\right(x-x_{\ast })/\mathit{\lambda })$ with a constant C₀ and a small ϵ≪1. Then the change in the steady-state solution u₁ at the surface is according to Eq. (10)

$$\begin{array}{}\text{(28)}& \mathit{\delta }{u}_{\mathrm{1}}(x_{\ast },\mathit{\lambda })=\underset{\mathrm{0}}{\overset{L}{\int }}\mathit{ϵ}{C}_{\mathrm{0}}\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{v}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }(x-x_{\ast })}{\mathit{\lambda }}}\right)\mathrm{d}x.\end{array}$$

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⁵. When λ is large compared to ϱ then

$$\begin{array}{ll}{\displaystyle }\mathit{\delta }{u}_{\mathrm{1}}(x_{\ast },\mathit{\lambda })\approx \mathit{\delta }{u}_{\mathrm{1},\mathrm{\infty }}\left(x_{\ast }\right)& {\displaystyle }=\underset{\mathit{\lambda }\to \mathrm{\infty }}{lim}\mathit{\delta }{u}_{\mathrm{1}}(x_{\ast },\mathit{\lambda })\\ \text{(29)}& {\displaystyle }& {\displaystyle }=\mathit{ϵ}{C}_{\mathrm{0}}\underset{\mathrm{0}}{\overset{L}{\int }}\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{v}\mathrm{d}x,\end{array}$$

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₁ and h is small. In Fig. 3, $\mathit{\delta }{u}_{\mathrm{1}}(x_{\ast },\mathit{\lambda })$ and $\mathit{\delta }{u}_{\mathrm{1},\mathrm{\infty }}\left(x_{\ast }\right)$ are compared at $x_{\ast }=\mathrm{0.9}\times {\mathrm{10}}^{\mathrm{6}}$. When λ>ϱ then $\mathit{\delta }{u}_{\mathrm{1}}(x_{\ast },\mathit{\lambda })\approx \mathit{\delta }{u}_{\mathrm{1},\mathrm{\infty }}\left(x_{\ast }\right)$. Suppose that $\mathit{\lambda }=\mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$. Then $\mathit{\delta }{u}_{\mathrm{1}}(x_{\ast },\mathit{\lambda })$ is about $\mathrm{0.02}\mathit{\delta }{u}_{\mathrm{1},\mathrm{\infty }}\left(x_{\ast }\right)$ and probably hard to observe and $\mathit{\delta }h(x_{\ast },\mathit{\lambda })\approx \mathrm{0.2}\mathit{\delta }{h}_{\mathrm{\infty }}\left(x_{\ast }\right)$. Similar conclusions are drawn theoretically in Gudmundsson (2003, 2008) using Fourier analysis and experimentally in Sun et al. (2014).

Figure 3 The response at Γ_s with different wavelengths λ in the perturbation of C in Eq. (28). (a) $\mathit{\delta }{u}_{\mathrm{1}}(x_{\ast },\mathit{\lambda })/\mathit{\delta }{u}_{\mathrm{1},\mathrm{\infty }}\left(x_{\ast }\right)$; (b) $\mathit{\delta }h(x_{\ast },\mathit{\lambda })/\mathit{\delta }{h}_{\mathrm{\infty }}\left(x_{\ast }\right)$.

We perform a pair of experiments to compare the results from perturbing the forward equation and the prediction by the adjoint solutions. A relative 1 % perturbation δC(x) is added at $x\in [\mathrm{0.9},\mathrm{1.0}]\times {\mathrm{10}}^{\mathrm{6}}$ m to the friction coefficient C(x). The differences between the forward FS solutions with and without the perturbation after 1 year are shown in Fig. 4 marked as “perturbed”. The predicted perturbations are computed from the solutions of the adjoint equation by varying x_∗ along the x axis and inserting it into Eq. (10). Each red dot in Fig. 4 corresponds to one single observation at x_∗. Both the u₁ and h predictions are in good agreement with the forward perturbations.

Figure 4 The changes on the horizontal velocity u₁ (a) and surface elevation h (b) in FS after 1 year with 1 % perturbation on C(x) at $x\in [\mathrm{0.9},\mathrm{1.0}]\times {\mathrm{10}}^{\mathrm{6}}$ m. Solid lines are the differences between the steady-state and perturbed transient solutions in Eq. (5). Red dots are the estimated perturbation using Eq. (10)

3.2 SSA

The same MISMIP benchmark experiment as in Sect. 3.1 is solved by the SSA on a one-dimensional uniform grid with mesh size Δx=1 km using standard finite difference methods implemented in MATLAB. The time derivatives are discretized by the implicit Euler method with a constant time step Δt=1 year as in Sect. 3.1. An upwind scheme is used for the spatial derivatives in the forward and adjoint advection equations to stabilize the numerical solutions. Replacing the Dirac delta with a Gaussian distribution function a few grid points wide in order to smoothen the observation function and avoid numerical oscillations in the solution has no major effect on the solutions.

The numerical solution of the forward SSA equations Eq. (17) is compared to the analytical approximations in the Appendix Eq. (A1) in Fig. 5. The detailed derivations of the analytical solutions in the Appendix are found in Cheng and Lötstedt (2019). 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 for the grounded ice. 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.

Figure 5 Comparison of the steady-state numerical solutions of the SSA velocity u and the thickness H in Eq. (17) (orange) and the analytical solutions in Eq. (A1) (blue).

The weight functions w_uC and w_hC in Fig. 6 have the same nonzero pattern as v since they are equal to −vu^m in Eq. (20). Each one of these weights w_uC or w_hC corresponds to the sensitivity of the observation at x_∗ with respect to the change in C(x), which is one row in the weight matrices W_uC or W_hC in Eqs. (23) and (24). The analytical weight functions in Eqs. (A3) and (A5) at $x_{\ast }=\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$ m are included in the steady-state figures for comparison. In the transient SSA simulations, the sensitivity is similar to those in the adjoint FS solutions in Fig. 2 increasing towards the grounding line. Such an increased sensitivity is also noted in Kyrke-Smith et al. (2018) and Leguy et al. (2014). However, in the steady-state cases, the weight functions indicate only an upstream effect of C(x). In other words, the perturbation in C(x) at point x can only influence the steady-state solutions to the left of this point. This is true as long as the effect of the grounding line migration is neglected. The δC weights for u responses are all negative, implying that an increase in C leads to decrease in u, but the steady-state surface elevation h rises when C is increased. The weights for the transient problem have a similar shape for the FS and SSA models in Figs. 2 and 6.

Figure 6 Comparison of the weights w_uC and w_hC for SSA in Eq. (19) for perturbations δC with m=1 at different observation points $x_{\ast }=\mathrm{0.25}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.5}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$, and 0.9×10⁶ (blue, orange, green, and pink). The dashed black lines in (c, d) are w_uC and w_hC computed from the analytical solutions of u in Eq. (A1) and v in Eqs. (A2) and (A4) at $x_{\ast }=\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$. (a, b) Transient simulations; (c, d) steady states. (a, c) w_uC for pointwise u response; (b, d) w_hC for pointwise h response.

The weight functions w_ub and w_hb for δb are localized at the observation position x_∗ in all the four cases in Fig. 7, which implies that the inverse problems may be well-posed. The dashed black lines in the two lower panels are the analytical expressions of the weight functions at $x_{\ast }=\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$ m in Eqs. (A3) and (A5) with a hat function of width 2Δx at the base to approximate the Dirac delta. The analytical solutions almost coincide with the numerical solutions. The steady-state weight functions are nonzero to the right of x_∗, corresponding to the integral in Eq. (A5). There is a detailed view of the steady-state δb weights for $x>x_{\ast }$ in Fig. 8. The weights of δb have similar structures to the δC weights. The analytical solutions in Eqs. (A3) and (A5) suggest that $w_{ub}/w_{uC}\approx w_{hb}/w_{hC}\approx (m+\mathrm{1})C/H$ for $x\ne x_{\ast }$.

Figure 7 Comparison of the weights w_ub and w_hb for SSA in Eq. (19) for perturbations δb at different observation points $x_{\ast }=\mathrm{0.25}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.5}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$, and 0.9×10⁶ (blue, orange, green, and pink). The dashed black lines in (c, d) are the weights of δb in Eqs. (A3) and (A5) at $x_{\ast }=\mathrm{0.7}\times {\mathrm{10}}^{\mathrm{6}}$. (a, b) Transient simulations; (c, d) steady states. (a, c) w_ub for pointwise u response; (b, d) w_hb for pointwise h response.

Figure 8 A close-up view of the steady-state weights in the lower panels of Fig. 7.

The same perturbation on C(x) as in Fig. 4 is imposed in the SSA simulations. The perturbed solutions after 1 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 1-year perturbation experiment in Fig. 10, the transient weight functions in the upper panels in Fig. 6 are used for the sensitivity estimates. The weight functions in the upper panels of Fig. 7 predict the response in Fig. 11.

The corresponding comparisons for the steady-state problem are made in Figs. 12 and 13 with the weights in the lower panels of Figs. 6 and 7. The analytical solutions of the steady-state perturbations from Eqs. (A3) and (A5) are shown with dashed black lines in these two figures.

The rapid change in δh in Figs. 10 and 11 is explained by the shape of the weight functions in the upper-right panels of Figs. 6 and 7. The weights can be approximated by $-\mathit{\theta }(x,t){\mathit{\delta }}^{\prime }(x-x_{\ast })$ for some θ>0. Then the surface response will be

$$\begin{array}{ll}{\displaystyle }\mathit{\delta }h\left(x_{\ast }\right)& {\displaystyle }=\underset{\mathrm{0}}{\overset{T}{\int }}\underset{\mathrm{0}}{\overset{L}{\int }}-\mathit{\theta }(x,t){\mathit{\delta }}^{\prime }(x-x_{\ast })\mathit{\delta }C\left(x\right)\mathrm{d}x\mathrm{d}t\\ {\displaystyle }& {\displaystyle }=\underset{\mathrm{0}}{\overset{T}{\int }}\left(\mathit{\theta }\mathit{\delta }C{)}^{\prime }\right(x_{\ast },t)\mathrm{d}t,\end{array}$$

where δC jumps discontinuously at $x=\mathrm{0.9}\times {\mathrm{10}}^{\mathrm{6}}$ and $x=\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{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 Eqs. (A3) and (A5). The discontinuities in the upper panel of Fig. 13 are caused by the jumps in δb at 0.9×10⁶ and 1.0×10⁶ and the first term in Eq. (A3). The jumps in δh in the lower panel of Fig. 13 are due to the first term in Eq. (A5).

All the predicted solutions from the adjoint SSA are in good agreement with the forward perturbation.

The inverse problem of the steady state for the friction coefficient may not be well-posed since the weights are all positive from x_∗ to x_GL. This is verified by checking the singular values of the sensitivity matrices W_uC and W_hC in Fig. 9, where the largest and smallest singular values of Σ_uC are 10⁻⁴ and 10⁻¹² with a large quotient σ_uC1∕σ_uCN and the span of the singular values of Σ_hC is from 10⁻⁴ to 10⁻⁸ (which is better).

The singular values of the sensitivity matrices W_ub and W_hb in Fig. 9 are in the interval 10⁻⁴ to 10⁻⁷ from large to small. They are better conditioned than the sensitivity matrices for C. In particular, Σ_hb (in pink) 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.

Figure 9 The singular values of the transfer matrices W_uC, W_hC, W_ub, and W_hb.

Good approximations of the sensitivity matrices W_uC, W_hC, W_ub, and W_hb are found in Eqs. (A3) and (A5) at given x_∗i and x_j as in Eq. (23). If the basal topography is unperturbed at the same x coordinate as the observation point such that $\mathit{\delta }{b}_{\ast }=\mathrm{0}$, then the contributions of δb and δC cannot be separated since they are both multiplied by the same weight except for a different scaling factor. This is in agreement with numerical investigations in Kyrke-Smith et al. (2018). It is shown in Cheng and Lötstedt (2019) that the perturbation in δu is proportional to the wavelength of δC. Perturbations with short wavelengths will not reach the surface. These conclusions are also drawn in numerical solutions of FS in Kyrke-Smith et al. (2018) and with transfer functions in the frequency space in Gudmundsson (2008). The perturbation in u due to δC increases with increasing u and decreasing H. The sensitivity of δu and δh behaves in a better way if the observation at x_∗ is above the perturbation at x in the topography in Eqs. (A3) and (A5). Then δb and its derivative directly affect the perturbations at the top of the ice. This is in agreement with the computed singular values in Fig. 9. This property is utilized in van Pelt et al. (2013) when the bottom topography is inferred from height data. Inferring the geometry of the base from such data is easier than inferring the slipperiness and C because of the first term in Eq. (A5) and w_hb in the right column of Fig. 7.

Figure 10 The changes in the horizontal velocity u (a) and surface elevation h (b) for SSA after 1 year with 1 % perturbation of C(x) in $x\in [\mathrm{0.9},\mathrm{1.0}]\times {\mathrm{10}}^{\mathrm{6}}$ m. Solid lines are the differences between the steady-state and the perturbed solutions in Eq. (13). Red dots represent the estimated perturbation using Eq. (15).

Figure 11 The changes in the horizontal velocity u (a) and surface elevation h (b) for SSA after 1 year with 0.01 m perturbation of b(x) in $x\in [\mathrm{0.9},\mathrm{1.0}]\times {\mathrm{10}}^{\mathrm{6}}$ m. Solid lines are the differences between the steady-state and the perturbed solutions in Eq. (13). Red dots represent the estimated perturbation using Eq. (15).

Figure 12 The changes in the horizontal velocity u (a) and surface elevation h (b) for SSA after 15 000 years (close to the steady state) with 1 % perturbation of C(x) in $x\in [\mathrm{0.9},\mathrm{1.0}]\times {\mathrm{10}}^{\mathrm{6}}$. Solid lines are the differences between the steady-state and perturbed solutions in Eq. (13). Red dots represent the estimated perturbation using Eq. (15).

Figure 13 The changes in the horizontal velocity u (a) and surface elevation h (b) for SSA after 15 000 years (close to the steady state) with 0.01 m perturbation of b(x) in $x\in [\mathrm{0.9},\mathrm{1.0}]\times {\mathrm{10}}^{\mathrm{6}}$. Solid lines are the differences between the steady-state and perturbed solutions in Eq. (13). Red dots represent the estimated perturbation using Eq. (15).

The solution of the adjoint equations is simplified in the comparison in Fig. 14. In MacAyeal (1993), two simplifications are made. Firstly, the adjoint viscosity $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ in Eq. (14) is approximated by the forward viscosity η in Eq. (11). The factor 1∕n in the viscosity in the 2-D stress equation Eq. (18) is then replaced by 1. Secondly, the thickness H is fixed and the advection equation for ψ is not solved, which is equivalent to ∇ψ=0 in the adjoint stress equation in Eq. (15). Perturbations are introduced in C and u is observed for the transient case as in Fig. 10. The perturbed forward solutions are compared to the predicted perturbations by the simplified adjoint SSA systems in Fig. 14, where the forward viscosity η is used in both cases. In the upper panel of Fig. 14, the two equations of ψ and v are solved. In the lower panel, the advection equation of ψ is excluded from the system. The differences are small in this case compared to the full adjoint solution used in Fig. 10. The reason is that $\mathit{\psi },{\mathit{\psi }}_x,$ and Hηu_x are small in Eq. (18).

Figure 14 The changes in the horizontal velocity u for SSA after 1 year with 1 % perturbation of C(x) in $x\in [\mathrm{0.9},\mathrm{1.0}]\times {\mathrm{10}}^{\mathrm{6}}$ m. Solid lines are the differences between the steady-state and the perturbed solutions in Eq. (13). Red dots represent the estimated perturbation using Eq. (15). (a) Forward viscosity. (b) Without advection equation.

The singular values of the transfer matrices corresponding to the two simplifications are displayed in Fig. 15, where the two transfer matrices are denoted by ${\stackrel{\mathrm{̃}}{\mathbf{W}}}_{uC}$ for the system coupling ψ and v and by ${\widehat{\mathbf{W}}}_{uC}$ for the adjoint equation without ψ with a fixed H. The singular values in ${\stackrel{\mathrm{̃}}{\mathbf{\Sigma }}}_{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 ${\widehat{\mathbf{W}}}_{uC}$ has a better conditioning than ${\stackrel{\mathrm{̃}}{\mathbf{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 for penalizing oscillatory behavior at the base as in Gagliardini et al. (2013) and Petra et al. (2012).

Figure 15 The singular values of the transfer matrices of SSA with simplifications from MacAyeal (1993). ${\stackrel{\mathrm{̃}}{\mathbf{\Sigma }}}_{uC}$ corresponds to the forward viscosity case and ${\widehat{\mathbf{\Sigma }}}_{uC}$ is from the adjoint SSA without coupling to the ψ equation.

4 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 $({\mathit{u}}^i,p^i,h^i),i=\mathrm{1},\mathrm{2},\mathrm{\dots },K,$ to be able to compute the adjoint solutions $({\mathit{v}}^i,q^i,{\mathit{\psi }}^i),i=K,K-\mathrm{1},\mathrm{\dots },\mathrm{1},$ which may be prohibitive in 3-D. Since the data to be stored in the SSA model are one dimension lower, we are able to solve the adjoint problem backward in time for any number of K. For a fair comparison, we show the results for one time step with SSA in this paper.

The equations for the adjoints of FS and SSA in Eqs. (8) and (15) are generally valid for an ice sheet in 3-D and have to be solved numerically. The problem with the storage of the forward solution is the same as in adaptive mesh refinement, where the time step and the mesh are adapted to satisfy bounds on the numerical error. Selected forward solutions in time are saved for the adjoint solution to reduce the storage requirements. Missing values are interpolated in time and the sensitivity integral in Eqs. (10) and (16) is computed successively when the adjoint solution is advanced backward in time.

The solutions of the horizontal velocity u and the height h with perturbations in C in the transient FS and SSA models are similar in Figs. 4 and 10. The weights in the upper panels in Figs. 2 and 6 are similar, too. The solutions to the forward equations are also close in the chosen MISMIP configuration. The reason is that the sliding on the ground in the FS model is considerable, making SSA a good approximation of FS.

There are many discussions regarding the choice of friction laws (see, e.g., Gladstone et al., 2017; Tsai et al., 2015; Brondex et al., 2017). However, assuming a spatial variability of the friction coefficient C(x) with a linear relation between the basal stress and velocity makes this numerical study independent of the friction law. The friction coefficient can be viewed as a linearization of the friction law, and a postprocessing procedure can retrieve the corresponding friction law.

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 appears 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 𝒲_F and the adjoint problem is 𝒲_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 ${\mathcal{W}}_A/{\mathcal{W}}_F\approx \mathrm{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 $(\mathrm{1}+M{\mathcal{W}}_A/{\mathcal{W}}_F)/(\mathrm{1}+N)$. Since k≥1, it is advantageous to choose the approach involving the adjoint if M<kN. Otherwise, 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 with a gradient method. Solving the adjoint equations is then always favorable.

5 Conclusions

The perturbations δu and δh in the velocity u and the height h at the ice surface are caused by perturbations δb and δC in the topography of the ice base b and the basal friction coefficient C. The sensitivities δu and δh to δb and δC are evaluated in 2-D by first solving the adjoint equations of the FS and SSA models including the advection equation for the height derived in Cheng and Lötstedt (2019). Then weight or transfer functions are determined for the relation between δu and δh at the surface and δb and δC at the base. The predictions of δu and δh with the weights are compared to explicit calculations of perturbed u and h at the surface with good agreement. It is shown in Cheng and Lötstedt (2019) that if the base perturbations are time dependent then it is necessary to have time-dependent weight functions to obtain the correct behavior at the top of the ice.

Both the height and the stress equations and their adjoints are solved to find the weight functions here. The inverse problem at steady state to infer C from observations of u is usually solved for a fixed ice geometry and with only the stress equation and its adjoint (see, e.g., MacAyeal, 1993; Petra et al., 2012). This is possible since the adjoint height ψ is small when the horizontal 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 proportional to the velocity u and 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. This is observed numerically in Durand et al. (2011).

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 $\mathit{\lambda }=\mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ m (about 10H) in C is damped by 0.2 in δh and 0.02 in δu compared to a wavelength λ>10⁵ where there is no damping due to λ. This is in agreement with the theory in Gudmundsson (2008).

The perturbations in u and h in the steady state of the SSA model consist of a direct effect from δb at the observation point and a nonlocal effect of δb and δC in Figs. 6 and 7. It follows from the analytical solution in Eq. (A3) that we cannot distinguish between the nonlocal contributions of δb and δC in the integral to δu. The same conclusion about the nonlocal perturbations holds for δh in Eq. (A5). This is also an observation in Kyrke-Smith et al. (2018).

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

$${\displaystyle }{\displaystyle }\mathrm{1}.\mathit{\delta }b\to \mathit{\delta }h,\mathrm{2}.\mathit{\delta }b\to \mathit{\delta }u,\mathrm{3}.\mathit{\delta }C\to \mathit{\delta }h,\mathrm{4}.\mathit{\delta }C\to \mathit{\delta }u.$$

In the past, the coupling between δu and δC is most frequently used for inference of C from velocity data, but adding height data would improve the robustness of the inference. The approximated analytical transfer functions for SSA, yielding explicit dependence of the parameters, have the same properties as above in which the observed velocity u and height h are more sensitive to perturbations δb than δC.

Appendix A: Some equations

Detailed derivations of the formulas are found in Cheng and Lötstedt (2019). A variable with index ∗ is evaluated at x_∗.

A1 The forward steady-state SSA solution

The analytical steady-state solution to the forward Eq. (17) without considering the viscosity terms is

$$\begin{array}{ll}{\displaystyle }H\left(x\right)& {\displaystyle }={\displaystyle }{\left(H_{\text{GL}}^{m+\mathrm{2}}+{\displaystyle \frac{m+\mathrm{2}}{m+\mathrm{1}}}{\displaystyle \frac{C{a}^m}{\mathit{\rho }g}}(x_{\text{GL}}^{m+\mathrm{1}}-x^{m+\mathrm{1}})\right)}^{\frac{\mathrm{1}}{m+\mathrm{2}}},\\ {\displaystyle }& {\displaystyle }\mathrm{0}\le x\le x_{\text{GL}},\\ {\displaystyle }H\left(x\right)& {\displaystyle }=H_{\text{GL}},x_{\text{GL}}<x<L,\\ \text{(A1)}& {\displaystyle }u\left(x\right)& {\displaystyle }={\displaystyle \frac{ax}{H}},\mathrm{0}\le x\le x_{\text{GL}},u\left(x\right)={\displaystyle \frac{ax}{H_{\text{GL}}}},x_{\text{GL}}<x<L,\end{array}$$

where H_GL is the thickness of the ice at the grounding line x_GL.

A2 The adjoint steady-state SSA solutions

The analytical steady-state solutions of the SSA adjoint Eq. (18) with observation of u at x_∗ are

$$\begin{array}{ll}{\displaystyle }\mathit{\psi }\left(x\right)& {\displaystyle }={\displaystyle \frac{C{a}^m{x}_{*}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{3}}}}\left(x_{\text{GL}}^m-x^m\right),x_{\ast }<x\le x_{\text{GL}},\\ {\displaystyle }\mathit{\psi }\left(x\right)& {\displaystyle }={\displaystyle }-{\displaystyle \frac{\mathrm{1}}{H_{\ast }}}+{\displaystyle \frac{C{a}^m{x}_{*}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{3}}}}\left(x_{\text{GL}}^m-x_{\ast }^m\right),\mathrm{0}\le x<x_{\ast },\\ {\displaystyle }v\left(x\right)& {\displaystyle }={\displaystyle \frac{a{x}_{*}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{3}}}}{H}^m,x_{\ast }<x\le x_{\text{GL}},\\ \text{(A2)}& {\displaystyle }v\left(x\right)& {\displaystyle }=\mathrm{0},\mathrm{0}\le x<x_{\ast },\end{array}$$

where H_∗ is the thickness of the ice at x_∗. The corresponding perturbation δu_∗ in Eq. (20) has the weights for δC and δb as

$$\begin{array}{ll}{\displaystyle }\mathit{\delta }{u}_{\ast }& {\displaystyle }={\displaystyle }\underset{\mathrm{0}}{\overset{x_{\text{GL}}}{\int }}({\mathit{\psi }}_{x}u+v_x\mathit{\eta }{u}_x+v\mathit{\rho }g{h}_x)\mathit{\delta }b-v{u}^m\mathit{\delta }C\mathrm{d}x\\ \text{(A3)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{u_{\ast }}{H_{\ast }}}\mathit{\delta }{b}_{\ast }-{\displaystyle \frac{u_{\ast }}{H_{\ast }}}\underset{x_{\ast }}{\overset{x_{\text{GL}}}{\int }}{\displaystyle \frac{C(ax{)}^m}{\mathit{\rho }g{H}_{*}^{m+\mathrm{1}}}}\left((m+\mathrm{1}){\displaystyle \frac{\mathit{\delta }b}{H}}+{\displaystyle \frac{\mathit{\delta }C}{C}}\right)\mathrm{d}x.\end{array}$$

If h is observed at x_∗, then

$$\begin{array}{ll}{\displaystyle }\mathit{\psi }\left(x\right)& {\displaystyle }={\displaystyle }-{\displaystyle \frac{C{a}^{m-\mathrm{1}}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{1}}}}\left(x_{\text{GL}}^m-x^m\right),x_{\ast }<x\le x_{\text{GL}},\\ {\displaystyle }\mathit{\psi }\left(x\right)& {\displaystyle }={\displaystyle }-{\displaystyle \frac{C{a}^{m-\mathrm{1}}}{\mathit{\rho }g{H}_{*}^{m+\mathrm{1}}}}\left(x_{\text{GL}}^m-x_{\ast }^m\right)-{\displaystyle \frac{\mathit{\delta }(x-x_{\ast }){\mathit{\eta }}_{\ast }}{n\mathit{\rho }g{H}_{\ast }}},\mathrm{0}\le x\le x_{\ast },\\ {\displaystyle }v\left(x\right)& {\displaystyle }={\displaystyle }-{\displaystyle \frac{H^m}{\mathit{\rho }g{H}_{*}^{m+\mathrm{1}}}},x_{\ast }<x\le x_{\text{GL}},\\ \text{(A4)}& {\displaystyle }v\left(x\right)& {\displaystyle }=\mathrm{0},\mathrm{0}\le x<x_{\ast }.\end{array}$$

The weights for δC and δb in Eq. (19) for the perturbation on h_∗ are

$$\begin{array}{ll}{\displaystyle }\mathit{\delta }{h}_{\ast }=& {\displaystyle \frac{{\mathit{\eta }}_{\ast }}{n\mathit{\rho }g{H}_{\ast }}}\left(u\mathit{\delta }b{)}_x\right(x_{\ast })\\ \text{(A5)}& {\displaystyle }& {\displaystyle }+{\displaystyle }\underset{x_{\ast }}{\overset{x_{\text{GL}}}{\int }}{\displaystyle \frac{C(ax{)}^m}{\mathit{\rho }g{H}_{*}^{m+\mathrm{1}}}}\left((m+\mathrm{1}){\displaystyle \frac{\mathit{\delta }b}{H}}+{\displaystyle \frac{\mathit{\delta }C}{C}}\right)\mathrm{d}x.\end{array}$$

A3 The singular value decomposition (SVD)

The SVD factorizes a matrix A in the following way (see Golub and Van Loan, 1989):

$$\begin{array}{}\text{(A6)}& \mathbf{A}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T.\end{array}$$

If A is an M×N matrix then U is an M×M matrix, Σ an M×N matrix, and V an N×N matrix. The singular values σ_i are nonnegative and ordered from large to small for increasing i and $i=\mathrm{1},\mathrm{2},\mathrm{\dots },min(M,N)$. They form the diagonal of the diagonal matrix Σ with Σ_ii=σ_i. The other two matrices are orthogonal, satisfying U^TU=I and V^TV=I. The generalized inverse Σ⁻¹ of Σ is an N×M matrix with ${\mathit{\sigma }}_i^{-\mathrm{1}}$ (if σ_i is positive) on the diagonal and 0 elsewhere.