- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer review
- For authors
- For reviewers
- EGU publications

Journal cover
Journal topic
**The Cryosphere**
An interactive open-access journal of the European Geosciences Union

Journal topic

- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer review
- For authors
- For reviewers
- EGU publications

**Research article**
17 Feb 2020

**Research article** | 17 Feb 2020

Parameter sensitivity analysis of dynamic ice sheet models – numerical computations

- Department of Information Technology, Uppsala University, P.O. Box 337, 75105 Uppsala, Sweden

- Department of Information Technology, Uppsala University, P.O. Box 337, 75105 Uppsala, Sweden

**Correspondence**: Gong Cheng (cheng.gong@it.uu.se)

**Correspondence**: Gong Cheng (cheng.gong@it.uu.se)

Abstract

Back to toptop
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.

Download & links

How to cite

Back to top
top
How to cite.

Cheng, G. and Lötstedt, P.: Parameter sensitivity analysis of dynamic ice sheet models – numerical computations, The Cryosphere, 14, 673–691, https://doi.org/10.5194/tc-14-673-2020, 2020.

1 Introduction

Back to toptop
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

$$\begin{array}{ll}{\displaystyle}& {\displaystyle}(\mathit{a}\otimes \mathit{c}{)}_{ij}={a}_{i}{c}_{j},\phantom{\rule{1em}{0ex}}\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},\phantom{\rule{1em}{0ex}}\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

Back to toptop
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).

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**,

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 Γ

$$\begin{array}{ll}{\displaystyle}\mathrm{\Omega}& {\displaystyle}=\mathit{\left\{}\mathit{x}\right|(x,y)\in \mathit{\omega},\phantom{\rule{0.125em}{0ex}}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},\phantom{\rule{0.125em}{0ex}}z=h(x,y,t)\mathit{\}},\\ {\displaystyle}{\mathrm{\Gamma}}_{b}& {\displaystyle}=\mathit{\left\{}\mathit{x}\right|(x,y)\in \mathit{\omega},\phantom{\rule{0.125em}{0ex}}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},\phantom{\rule{0.125em}{0ex}}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},\phantom{\rule{0.125em}{0ex}}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},\phantom{\rule{0.125em}{0ex}}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},\phantom{\rule{0.25em}{0ex}}\mathit{\omega}=[\mathrm{0},L],\phantom{\rule{0.25em}{0ex}}{\mathit{\gamma}}_{u}=\mathrm{0},$ and *γ*_{d}=*L*.

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}),\phantom{\rule{0.25em}{0ex}}\mathit{\eta}\left(\mathit{u}\right)=\frac{\mathrm{1}}{\mathrm{2}}{A}^{-{\scriptscriptstyle \frac{\mathrm{1}}{n}}}\left(\text{tr}{\mathbf{D}}^{\mathrm{2}}\right(\mathit{u}){)}^{\mathit{\nu}},\phantom{\rule{0.25em}{0ex}}\mathit{\nu}=\frac{\mathrm{1}-n}{\mathrm{2}n}.\end{array}\end{array}$$

The trace of **D**^{2} is tr**D**^{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

$$\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

$$\begin{array}{ll}{\displaystyle}& {\displaystyle}{h}_{t}+\mathit{h}\cdot \mathit{u}=a,\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Gamma}}_{s},\\ {\displaystyle}& {\displaystyle}h(\mathit{x},\mathrm{0})={h}_{\mathrm{0}}\left(\mathit{x}\right),\phantom{\rule{0.25em}{0ex}}\mathit{x}\in \mathit{\omega},\phantom{\rule{1em}{0ex}}h(\mathit{x},t)={h}_{\mathit{\gamma}}(\mathit{x},t),\phantom{\rule{0.25em}{0ex}}\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},\phantom{\rule{0.25em}{0ex}}\text{in}\phantom{\rule{0.25em}{0ex}}\mathrm{\Omega}\left(t\right),\phantom{\rule{1em}{0ex}}\mathit{\sigma}\mathit{n}=\mathbf{0},\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Gamma}}_{s},\\ \text{(5)}& {\displaystyle}& {\displaystyle}\mathbf{T}\mathit{\sigma}\mathit{n}=-Cf\left(\mathbf{T}\mathit{u}\right)\mathbf{T}\mathit{u},\phantom{\rule{1em}{0ex}}\mathit{n}\cdot \mathit{u}=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Gamma}}_{b}.\end{array}$$

The initial data for *h* are *h*_{0}(** x**), and

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

The Dirichlet boundary conditions of ** u** on Γ

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)\phantom{\rule{0.125em}{0ex}}\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

$$}{\displaystyle}\mathcal{F}=\underset{{\mathrm{\Gamma}}_{s}}{\int}F\left(\mathit{u}\right)\phantom{\rule{0.125em}{0ex}}\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}\mathcal{F}=\underset{{\mathrm{\Gamma}}_{s}}{\int}F\left(h\right)\phantom{\rule{0.125em}{0ex}}\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**,

The adjoint FS equations form a system of PDEs for the adjoint height *ψ*, the adjoint velocity ** v**, and the adjoint pressure

$$\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},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Gamma}}_{s},\\ {\displaystyle}& {\displaystyle}\mathit{\psi}(\mathit{x},T)=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\mathit{\psi}(\mathit{x},t)=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Gamma}}_{d},\\ {\displaystyle}& {\displaystyle}-\mathrm{\nabla}\cdot \stackrel{\mathrm{\u0303}}{\mathit{\sigma}}(\mathit{v},q)=-\mathrm{\nabla}\cdot \left(\mathrm{2}\stackrel{\mathrm{\u0303}}{\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},\phantom{\rule{0.25em}{0ex}}\text{in}\phantom{\rule{0.25em}{0ex}}\mathrm{\Omega}\left(t\right),\\ {\displaystyle}& {\displaystyle}\stackrel{\mathrm{\u0303}}{\mathit{\sigma}}(\mathit{v},q)\mathit{n}=-({F}_{\mathit{u}}+\mathit{\psi}\mathit{h}),\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Gamma}}_{s},\\ {\displaystyle}& {\displaystyle}\mathbf{T}\stackrel{\mathrm{\u0303}}{\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},\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Gamma}}_{b},\\ \text{(8)}& {\displaystyle}& {\displaystyle}\mathit{n}\cdot \mathit{v}=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\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{\u0303}}{\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{\u0303}}{\mathit{\sigma}}(\mathit{v},q)& {\displaystyle}=\mathrm{2}\stackrel{\mathrm{\u0303}}{\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 *i**j**k**l* 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}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathit{\delta}C\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{x}\phantom{\rule{0.125em}{0ex}}\mathrm{d}t,\end{array}$$

involving the tangential projections of the forward and adjoint velocities **T**** u** and

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.

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.

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}\phantom{\rule{0.125em}{0ex}}\mathbf{I}$.
The stress tensor ** ς**(

$$\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

The ice dynamics system is

$$\begin{array}{ll}{\displaystyle}& {\displaystyle}{h}_{t}+\mathrm{\nabla}\cdot \left(\mathit{u}H\right)=a,\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathrm{0}\le t\le T,\phantom{\rule{0.25em}{0ex}}\mathit{x}\in \mathit{\omega},\\ {\displaystyle}& {\displaystyle}h(\mathit{x},\mathrm{0})={h}_{\mathrm{0}}\left(\mathit{x}\right),\phantom{\rule{0.25em}{0ex}}\mathit{x}\in \mathit{\omega},\phantom{\rule{1em}{0ex}}h(\mathit{x},t)={h}_{\mathit{\gamma}}(\mathit{x},t),\phantom{\rule{0.125em}{0ex}}\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,\phantom{\rule{0.25em}{0ex}}\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},\phantom{\rule{0.25em}{0ex}}\mathit{x}\in {\mathit{\gamma}}_{g},\\ \text{(13)}& {\displaystyle}& {\displaystyle}\mathit{t}\cdot \mathit{\varsigma}\mathit{n}=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\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**⋅

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

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{\u0303}}{\mathit{\eta}}$ and adjoint stress $\stackrel{\mathrm{\u0303}}{\mathit{\varsigma}}$ are defined by

$$\begin{array}{ll}{\displaystyle}\stackrel{\mathrm{\u0303}}{\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{\u0303}}{\mathit{\varsigma}}\left(\mathit{v}\right)& {\displaystyle}=\mathrm{2}H\stackrel{\mathrm{\u0303}}{\mathit{\eta}}\left(\mathit{u}\right)\star \mathbf{B}\left(\mathit{v}\right);\end{array}$$

cf. $\stackrel{\mathrm{\u0303}}{\mathit{\eta}}$ and $\stackrel{\mathrm{\u0303}}{\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},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{in}\phantom{\rule{0.25em}{0ex}}\mathit{\omega},\\ {\displaystyle}& {\displaystyle}\mathit{\psi}(\mathit{x},T)=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{in}\phantom{\rule{0.25em}{0ex}}\mathit{\omega},\phantom{\rule{1em}{0ex}}\mathit{\psi}(\mathit{x},t)=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathit{\gamma}}_{w},\\ {\displaystyle}& {\displaystyle}\mathrm{\nabla}\cdot \stackrel{\mathrm{\u0303}}{\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},\phantom{\rule{1em}{0ex}}\text{in}\phantom{\rule{0.25em}{0ex}}\mathit{\omega},\\ {\displaystyle}& {\displaystyle}\mathit{t}\cdot \stackrel{\mathrm{\u0303}}{\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},\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathit{\gamma}}_{g},\\ {\displaystyle}& {\displaystyle}\mathit{t}\cdot \stackrel{\mathrm{\u0303}}{\mathit{\varsigma}}\left(\mathit{v}\right)\mathit{n}=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}{\mathit{\gamma}}_{w},\\ \text{(15)}& {\displaystyle}& {\displaystyle}\mathit{n}\cdot \mathit{v}=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\text{on}\phantom{\rule{0.25em}{0ex}}\mathit{\gamma}.\end{array}$$

Compared to Eq. (8), the advection equation depends on ** v**, and the influence of

$$}{\displaystyle}{\mathbf{F}}_{\mathit{\omega}}\left(\mathit{u}\right)={\displaystyle \frac{m-\mathrm{1}}{\mathit{u}\cdot \mathit{u}}}\mathit{u}\otimes \mathit{u},\phantom{\rule{1em}{0ex}}{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})\phantom{\rule{0.125em}{0ex}}\mathit{\delta}b\\ {\displaystyle}& {\displaystyle}-f\left(\mathit{u}\right)\mathit{u}\cdot \mathit{v}\phantom{\rule{0.125em}{0ex}}\mathit{\delta}C\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{x}\phantom{\rule{0.125em}{0ex}}\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}\phantom{\rule{0.125em}{0ex}}\mathit{t}\cdot \mathit{v}\phantom{\rule{0.125em}{0ex}}\mathit{\delta}{C}_{\mathit{\gamma}}\phantom{\rule{0.125em}{0ex}}\mathrm{d}s\phantom{\rule{0.125em}{0ex}}\mathrm{d}t.\end{array}$$

In the 2-D model, *u*_{2}=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*_{1} be independent of *y* and taking *u*_{2}=0.
The notation is simplified if we let *u*=*u*_{1} and *v*=*v*_{1}.
The forward equations follow from Eq. (13):

$$\begin{array}{ll}{\displaystyle}& {\displaystyle}{h}_{t}+(uH{)}_{x}=a,\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathrm{0}\le t\le T,\phantom{\rule{0.25em}{0ex}}\mathrm{0}\le x\le L,\\ {\displaystyle}& {\displaystyle}h(x,\mathrm{0})={h}_{\mathrm{0}}\left(x\right),\phantom{\rule{0.25em}{0ex}}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},\phantom{\rule{0.25em}{0ex}}\mathrm{0}\le x\le L,\\ \text{(17)}& {\displaystyle}& {\displaystyle}u(\mathrm{0},t)={u}_{L}\left(t\right),\phantom{\rule{0.25em}{0ex}}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 *C**f*(*u*)*u*=*C**u*^{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,\phantom{\rule{0.25em}{0ex}}\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},\phantom{\rule{0.25em}{0ex}}\mathit{\psi}(L,t)=\mathrm{0},\phantom{\rule{0.25em}{0ex}}v(\mathrm{0},t)=\mathrm{0},\phantom{\rule{0.25em}{0ex}}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})\phantom{\rule{0.125em}{0ex}}\mathit{\delta}b-vf\left(u\right)u\phantom{\rule{0.125em}{0ex}}\mathit{\delta}C\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathrm{d}x\phantom{\rule{0.125em}{0ex}}\mathrm{d}t.\end{array}$$

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\phantom{\rule{0.125em}{0ex}}\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},\phantom{\rule{0.25em}{0ex}}{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\phantom{\rule{0.125em}{0ex}}\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},\phantom{\rule{0.125em}{0ex}}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},\phantom{\rule{0.125em}{0ex}}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))\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}x,\\ \text{(22)}& {\displaystyle}{\mathit{\mu}}_{\mathrm{1}}=& {\displaystyle}\phantom{\rule{0.125em}{0ex}}\mathrm{0.5},\phantom{\rule{0.125em}{0ex}}{\mathit{\mu}}_{j}=\mathrm{1},j=\mathrm{2},\mathrm{3},\mathrm{\dots},N-\mathrm{1},\phantom{\rule{0.125em}{0ex}}{\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}),\phantom{\rule{0.25em}{0ex}}{W}_{uCij}={\mathit{\mu}}_{j}{w}_{uC}({x}_{\ast i},{x}_{j}),\\ {\displaystyle}& {\displaystyle}i=\mathrm{1},\mathrm{2},\mathrm{\dots},M,\phantom{\rule{0.25em}{0ex}}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}{\mathbf{W}}_{ub}={\mathbf{U}}_{ub}{\mathbf{\Sigma}}_{ub}{\mathbf{V}}_{ub}^{T},\phantom{\rule{0.25em}{0ex}}{\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

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.

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

$$\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

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

Back to toptop
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}$$

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}}\phantom{\rule{0.125em}{0ex}}{\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}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}/\mathrm{3}}\phantom{\rule{0.125em}{0ex}}\mathrm{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 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 *C**f*(** u**)=

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 Γ

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*_{1} 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 Γ

The adjoint solutions *v*_{1} at Γ_{b} of all the four cases are concentrated at the observation points.
The vertical component *v*_{2} shares the same feature as *v*_{1} due to the boundary condition $\mathit{n}\cdot \mathit{v}=\mathrm{0}$ on Γ_{b}.
Therefore, the weights **T**** u**⋅

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{\u03f5}{C}_{\mathrm{0}}\mathrm{cos}\left(\mathrm{2}\mathit{\pi}\right(x-{x}_{\ast})/\mathit{\lambda})$ 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)

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

The same relation holds for *δ**h*(*x*_{∗}) but with a different ** v**. Let

$$\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{\u03f5}{C}_{\mathrm{0}}\underset{\mathrm{0}}{\overset{L}{\int}}\mathbf{T}\mathit{u}\cdot \mathbf{T}\mathit{v}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\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*_{1} 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).

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*_{1} and *h* predictions are in good agreement with the forward perturbations.

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

The weight functions *w*_{uC} and *w*_{hC} in Fig. 6 have the same nonzero pattern as *v* since they are equal to −*v**u*^{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.

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}$.

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)\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\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)\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\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^{6} and 1.0×10^{6} 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^{−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) 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*.

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.

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{\u0303}}{\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).

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{\u0303}}{\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{\u0303}}{\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{\u0303}}{\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).

4 Discussion

Back to toptop
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}),\phantom{\rule{0.25em}{0ex}}i=\mathrm{1},\mathrm{2},\mathrm{\dots},K,$ to be able to compute the adjoint solutions $({\mathit{v}}^{i},{q}^{i},{\mathit{\psi}}^{i}),\phantom{\rule{0.25em}{0ex}}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

5 Conclusions

Back to toptop
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

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

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 10*H*) 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 *λ*. 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

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}\mathrm{1}.\phantom{\rule{0.25em}{0ex}}\mathit{\delta}b\to \mathit{\delta}h,\phantom{\rule{0.25em}{0ex}}\mathrm{2}.\phantom{\rule{0.25em}{0ex}}\mathit{\delta}b\to \mathit{\delta}u,\phantom{\rule{0.25em}{0ex}}\mathrm{3}.\phantom{\rule{0.25em}{0ex}}\mathit{\delta}C\to \mathit{\delta}h,\phantom{\rule{0.25em}{0ex}}\mathrm{4}.\phantom{\rule{0.25em}{0ex}}\mathit{\delta}C\to \mathit{\delta}u.$$

In the past, the coupling between *δ*** u** and

Appendix A: Some equations

Back to toptop
Detailed derivations of the formulas are found in Cheng and Lötstedt (2019).
A variable with index ∗ is evaluated at *x*_{∗}.

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}},\phantom{\rule{0.25em}{0ex}}{x}_{\text{GL}}<x<L,\\ \text{(A1)}& {\displaystyle}u\left(x\right)& {\displaystyle}={\displaystyle \frac{ax}{H}},\phantom{\rule{0.25em}{0ex}}\mathrm{0}\le x\le {x}_{\text{GL}},\phantom{\rule{1em}{0ex}}u\left(x\right)={\displaystyle \frac{ax}{{H}_{\text{GL}}}},\phantom{\rule{0.25em}{0ex}}{x}_{\text{GL}}<x<L,\end{array}$$

where *H*_{GL} is the thickness of the ice at the grounding line *x*_{GL}.

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),\phantom{\rule{0.25em}{0ex}}{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),\phantom{\rule{0.25em}{0ex}}\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},\phantom{\rule{0.25em}{0ex}}{x}_{\ast}<x\le {x}_{\text{GL}},\\ \text{(A2)}& {\displaystyle}v\left(x\right)& {\displaystyle}=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\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})\phantom{\rule{0.125em}{0ex}}\mathit{\delta}b-v{u}^{m}\phantom{\rule{0.125em}{0ex}}\mathit{\delta}C\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\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)\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\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),\phantom{\rule{0.25em}{0ex}}{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}}},\phantom{\rule{0.25em}{0ex}}\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}}}},\phantom{\rule{0.25em}{0ex}}{x}_{\ast}<x\le {x}_{\text{GL}},\\ \text{(A4)}& {\displaystyle}v\left(x\right)& {\displaystyle}=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\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)\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathrm{d}x.\end{array}$$

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**^{T}**U**=**I** and **V**^{T}**V**=**I**. The generalized inverse **Σ**^{−1} of **Σ** is an *N*×*M* matrix with ${\mathit{\sigma}}_{i}^{-\mathrm{1}}$ (if *σ*_{i} is positive) on the diagonal and 0 elsewhere.

Code availability

Back to toptop
Code availability.

The FS equations are solved using Elmer/Ice version 8.4 (rev. f6bfdc9) with the scripts at: https://doi.org/10.5281/zenodo.3611158 (Cheng, 2020a). The forward and adjoint SSA solvers are implemented in MATLAB. The code is available at: https://doi.org/10.5281/zenodo.3611154 (Cheng, 2020b).

Author contributions

Back to toptop
Author contributions.

GC contributed most of the computations, and GC and PL contributed equally to the theory and the writing of the paper.

Competing interests

Back to toptop
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements

Back to toptop
Acknowledgements.

Thomas Zwinger helped with the adjoint FS solver in Elmer/Ice. Comments by Lina von Sydow helped us improve a draft of the paper, and Murtazo Nazarov explained the algorithm for adaptive mesh refinement.

Financial support

Back to toptop
Financial support.

This research has been supported by a Svenska Forskningsrådet Formas grant (no. 2017-00665) to Nina Kirchner and the Swedish strategic research program eSSENCE.

Review statement

Back to toptop
Review statement.

This paper was edited by Alexander Robinson and reviewed by two anonymous referees.

References

Back to toptop
Baiocchi, C., Brezzi, F., and Franca, L. P.: Virtual bubbles and Galerkin-least-squares type methods (Ga. LS), Comp. Meth. Appl. Mech. Eng., 105, 125–141, 1993. a

Brondex, J., Gagliardini, O., Gillet-Chaulet, F., and Durand, G.: Sensitivity of grounding line dynamics to the choice of the friction law, J. Glaciol., 63, 854–866, 2017. a

Bulthuis, K., Arnst, M., Sun, S., and Pattyn, F.: Uncertainty quantification of the multi-centennial response of the Antarctic ice sheet to climate change, The Cryosphere, 13, 1349–1380, https://doi.org/10.5194/tc-13-1349-2019, 2019. a

Cheng, G.: Numerical experiments for FS adjoint, Zenodo, https://doi.org/10.5281/zenodo.3611158, 2020a. a

Cheng, G.: Numerical experiments for SSA adjoint, Zenodo, https://doi.org/10.5281/zenodo.3611154, 2020b. a

Cheng, G. and Lötstedt, P.: Parameter sensitivity analysis of dynamic ice sheet models, Tech. Rep., arXiv:1906.08197, 2019. a, b, c, d, e, f, g, h, i, j, k, l

Durand, G., Gagliardini, O., Favier, L., Zwinger, T., and Le Meur, E.: Impact of bedrock description on modeling ice sheet dynamics, Geophys. Res. Lett., 38, L20501, https://doi.org/10.1029/2011GL048892, 2011. a, b

Gagliardini, O., Zwinger, T., Gillet-Chaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a new-generation ice sheet model, Geosci. Model Dev., 6, 1299–1318, https://doi.org/10.5194/gmd-6-1299-2013, 2013. a, b, c, d

Gillet-Chaulet, F., Gagliardini, O., Seddik, H., Nodet, M., Durand, G., Ritz, C., Zwinger, T., Greve, R., and Vaughan, D. G.: Greenland ice sheet contribution to sea-level rise from a new-generation ice-sheet model, The Cryosphere, 6, 1561–1576, https://doi.org/10.5194/tc-6-1561-2012, 2012. a

Gillet-Chaulet, F., Durand, G., Gagliardini, O., Mosbeux, C., Mouginot, J., Rémy, F., and Ritz, C.: Assimilation of surface velocities acquired between 1996 and 2010 to constrain the form of the basal friction law under Pine Island Glacier, Geophys. Res. Lett., 43, 10311–10321, 2016. a, b

Gladstone, R. M., Warner, R. C., Galton-Fenzi, B. K., Gagliardini, O., Zwinger, T., and Greve, R.: Marine ice sheet model performance depends on basal sliding physics and sub-shelf melting, The Cryosphere, 11, 319–329, https://doi.org/10.5194/tc-11-319-2017, 2017. a

Goldberg, D., Heimbach, P., Joughin, I., and Smith, B.: Committed retreat of Smith, Pope, and Kohler Glaciers over the next 30 years inferred by transient model calibration, The Cryosphere, 9, 2429–2446, https://doi.org/10.5194/tc-9-2429-2015, 2015. a

Golub, G. H. and Van Loan, C. F.: Matrix Computations, Johns Hopkins University Press, Baltimore, 2nd edn., 1989. a, b, c

Greve, R. and Blatter, H.: Dynamics of Ice Sheets and Glaciers, Advances in
Geophysical and Environmental Mechanics and Mathematics (AGEM^{2}), Springer, Berlin, 2009. a, b, c, d

Gudmundsson, G. H.: Transmission of basal variability to glacier surface, J. Geophys. Res., 108, 2253, https://doi.org/10.1029/2002JB002107, 2003. a, b, c, d

Gudmundsson, G. H.: Analytical solutions for the surface response to small amplitude perturbations in boundary data in the shallow-ice-stream approximation, The Cryosphere, 2, 77–93, https://doi.org/10.5194/tc-2-77-2008, 2008. a, b, c, d, e, f

Isaac, T., Petra, N., Stadler, G., and Ghattas, O.: Scalable and efficient algorithms for the propagation of uncertainty from data through inference to prediction for large-scale problems with application to flow of the Antarctic ice sheet, J. Comput. Phys., 296, 348–368, 2015. a, b

Jay-Allemand, M., Gillet-Chaulet, F., Gagliardini, O., and Nodet, M.: Investigating changes in basal conditions of Variegated Glacier prior to and during its 1982–1983 surge, The Cryosphere, 5, 659–672, https://doi.org/10.5194/tc-5-659-2011, 2011. a

Kyrke-Smith, T. M., Gudmundsson, G. H., and Farrell, P. E.: Relevance of detail in basal topography for basal slipperiness inversions: a case study on Pine Island Glacier, Antarctica, Frontiers Earth Sci., 6, 33, https://doi.org/10.3389/feart.2018.00033, 2018. a, b, c, d, e

Larour, E., Utke, J., Csatho, B., Schenk, A., Seroussi, H., Morlighem, M., Rignot, E., Schlegel, N., and Khazendar, A.: Inferred basal friction and surface mass balance of the Northeast Greenland Ice Stream using data assimilation of ICESat (Ice Cloud and land Elevation Satellite) surface altimetry and ISSM (Ice Sheet System Model), The Cryosphere, 8, 2335–2351, https://doi.org/10.5194/tc-8-2335-2014, 2014. a

Leguy, G. R., Asay-Davis, X. S., and Lipscomb, W. H.: Parameterization of basal friction near grounding lines in a one-dimensional ice sheet model, The Cryosphere, 8, 1239–1259, https://doi.org/10.5194/tc-8-1239-2014, 2014. a

MacAyeal, D. R.: Large-scale ice flow over a viscous basal sediment: Theory and application to Ice Stream B, Antarctica, J. Geophys. Res., 94, 4071–4078, 1989. a, b

MacAyeal, D. R.: A tutorial on the use of control methods in ice sheet modeling, J. Glaciol., 39, 91–98, 1993. a, b, c, d

Martin, N. and Monnier, J.: Adjoint accuracy for the full Stokes ice flow model: limits to the transmission of basal friction variability to the surface, The Cryosphere, 8, 721–741, https://doi.org/10.5194/tc-8-721-2014, 2014. a

Monnier, J. and des Boscs, P.-E.: Inference of the bottom properties in shallow ice approximation models, Inverse Problems, 33, 115001, https://doi.org/10.1088/1361-6420/aa7b92, 2017. a

Morlighem, M., Seroussi, H., Larour, E., and Rignot, E.: Inversion of basal friction in Antarctica using exact and incomplete adjoints of a high-order model, J. Geophys. Res.-Earth Surf., 118, 1–8, 2013. a

Mosbeux, C., Gillet-Chaulet, F., and Gagliardini, O.: Comparison of adjoint and nudging methods to initialise ice sheet model basal conditions, Geosci. Model Dev., 9, 2549–2562, https://doi.org/10.5194/gmd-9-2549-2016, 2016. a

Pattyn, F., Schoof, C., Perichon, L., Hindmarsh, R. C. A., Bueler, E., de Fleurian, B., Durand, G., Gagliardini, O., Gladstone, R., Goldberg, D., Gudmundsson, G. H., Huybrechts, P., Lee, V., Nick, F. M., Payne, A. J., Pollard, D., Rybak, O., Saito, F., and Vieli, A.: Results of the Marine Ice Sheet Model Intercomparison Project, MISMIP, The Cryosphere, 6, 573–588, https://doi.org/10.5194/tc-6-573-2012, 2012. a, b

Petra, N., Zhu, H., Stadler, G., Hughes, T. J. R., and Ghattas, O.: An inexact Gauss-Newton method for inversion of basal sliding and rheology parameters in a nonlinear Stokes ice sheet model, J. Glaciol., 58, 889–903, 2012. a, b, c, d, e, f, g, h

Schannwell, C., Drews, R., Ehlers, T. A., Eisen, O., Mayer, C., and Gillet-Chaulet, F.: Kinematic response of ice-rise divides to changes in ocean and atmosphere forcing, The Cryosphere, 13, 2673–2691, https://doi.org/10.5194/tc-13-2673-2019, 2019. a

Schlegel, N.-J., Seroussi, H., Schodlok, M. P., Larour, E. Y., Boening, C., Limonadi, D., Watkins, M. M., Morlighem, M., and van den Broeke, M. R.: Exploration of Antarctic Ice Sheet 100-year contribution to sea level rise and associated model uncertainties using the ISSM framework, The Cryosphere, 12, 3511–3534, https://doi.org/10.5194/tc-12-3511-2018, 2018. a

Schoof, C.: Ice-sheet acceleration driven by melt supply variability, Nature, 468, 803–806, 2010. a, b

Seddik, H., Greve, R., Sakakibara, D., Tsutaki, S., Minowa, M., and Sugiyama, S.: Response of the flow dynamics of Bowdoin Glacier, northwestern Greenland, to basal lubrication and tidal forcing, J. Glaciol., 65, 225–238, https://doi.org/10.1017/jog.2018.106, 2019. a

Sergienko, O. and Hindmarsh, R. C. A.: Regular patterns in frictional resistance of ice-stream beds seen by surface data inversion, Science, 342, 1086–1089, 2013. a

Shannon, S. R., Payne, A. J., Bartholomew, I. D., van den Broeke, M. R., Edwards, T. L., Fettweis, X., Gagliardini, O., Gillet-Chaulet, F., Goelzer, H., Hoffman, M. J., Huybrechts, P., Mair, D. W. F., Nienow, P. W., Perego, M., Price, S. F., Smeets, C. J. P. P., Sole, A. J., van de Wal, R. S. W., and Zwinger, T.: Enhanced basal lubrication and the contribution of the Greenland ice sheet to future sea–level rise, P. Natl. Acad. Sci. USA, 110, 14156–14161, 2013. a, b

Smith, R. C.: Uncertainty Quantification. Theory, Implementation, and Applications, Society for Industrial and Applied Mathematics, Philadelphia, 2014. a, b

Sole, A. J., Mair, D. W. F., Nienow, P. W., Bartholomew, I. D., King, I. D., Burke, M. A., and Joughin, I.: Seasonal speedup of a Greenland marine-terminating outlet glacier forced by surface melt-induced changes in subglacial hydrology, J. Geophys. Res., 116, F03014, https://doi.org/10.1029/2010JF001948, 2011. a

Sun, S., Cornford, S. L., Liu, Y., and Moore, J. C.: Dynamic response of Antarctic ice shelves to bedrock uncertainty, The Cryosphere, 8, 1561–1576, https://doi.org/10.5194/tc-8-1561-2014, 2014. a

Tsai, V. C., Stewart, A. L., and Thompson, A. F.: Marine ice-sheet profiles and stability under Coulomb basal conditions, J. Glaciol., 61, 205–215, 2015. a

Vallot, D., Pettersson, R., Luckman, A., Benn, D. I., Zwinger, T., van Pelt, W. J. J., Kohler, J., Schäfer, M., Claremar, B., and Hulton, N. R. J.: Basal dynamics of Kronebreen, a fast-flowing tidewater glacier in Svalbard: non-local spatio-temporal response to water input, J. Glaciol., 11, 179–190, 2017. a

van Pelt, W. J. J., Oerlemans, J., Reijmer, C. H., Pettersson, R., Pohjola, V. A., Isaksson, E., and Divine, D.: An iterative inverse method to estimate basal topography and initialize ice flow models, The Cryosphere, 7, 987–1006, https://doi.org/10.5194/tc-7-987-2013, 2013. a, b

Weertman, J.: On the sliding of glaciers, J. Glaciol., 3, 33–38, 1957. a

Short summary

We present a time-dependent inverse method for ice sheet modeling. By investigating the sensitivity of the observations of the velocity and the height at the surface to the basal conditions of the ice, we show that if the basal parameters are time dependent, then time cannot be ignored in the inversion. By looking at the numerical features, we conclude that adding the height information of an ice sheet in the velocity inversion procedure could improve the robustness of the inference.

We present a time-dependent inverse method for ice sheet modeling. By investigating the...

The Cryosphere

An interactive open-access journal of the European Geosciences Union