Inverting ice surface elevation and velocity for bed topography and slipperiness beneath Thwaites Glacier

There is significant uncertainty over how ice sheets and glaciers will respond to rising global temperatures. Limited knowledge of the topography and rheology of the ice–bed interface is a key cause of this uncertainty as models show that small changes in the bed can have a large influence on predicted rates of ice loss. Most of our detailed knowledge of bed topography comes from airborne and ground-penetrating radar observations. However, these direct observations are not spaced closely enough to meet the requirements of ice-sheet models, so interpolation and inversion methods are used to fill in the gaps. Here we present the results of a new inversion of surface elevation and velocity data over Thwaites Glacier, West Antarctica, for bed topography and slipperiness (i.e. the degree of basal slip for a given level of drag). The inversion is based on a steady-state linear perturbation analysis of the shallow-ice-stream equations. The method works by identifying disturbances to surface flow which are caused by obstacles or sticky patches in the bed and can therefore be applied wherever the shallow-ice-stream equations hold and where surface data are available, even where the ice thickness is not well known. We assess the performance of the inversion for topography with the available radar data. Although the topographic output from the inversion is less successful where the bed slopes steeply, it compares well with radar data from the central trunk of the glacier for medium-wavelength features (5–50 km). This method could therefore be useful as an independent test of other interpolation methods such as mass conservation and kriging. We do not have data to allow us to assess the success of the slipperiness results from our inversions, but we provide maps that may guide future seismic data collection across Thwaites Glacier. The methods presented here show significant promise for using high-resolution satellite datasets, calibrated by sparser field datasets, to generate high-resolution bed topography products across the ice sheets and therefore contribute to reduced uncertainty in predictions of future sea-level rise.

1 Introduction

Predicting the rate at which marine sectors of the West Antarctic Ice Sheet will retreat and contribute to globally rising sea levels is of increasing importance due to persistent climate forcing across the region over the last few decades (Scambos et al., 2017; Turner et al., 2017). Ice-sheet modelling studies emphasise the role of bed topography and rheology in understanding future ice loss (Durand et al., 2011; Parizek et al., 2013; Sun et al., 2014; Nias et al., 2016, 2018; Kyrke-Smith et al., 2018; Yu et al., 2018; Koellner et al., 2019). Bed topography is particularly important for marine-terminating glaciers, such as Thwaites Glacier in West Antarctica, which are vulnerable to marine ice-sheet instability (Weertman, 1974; Hughes, 1981; Schoof, 2007; Goldberg et al., 2009; Gudmundsson, 2013). However, bed topography constrained by geophysical surveying at the resolutions required for ice-sheet modelling (Durand et al., 2011; McCormack et al., 2018) is rarely available, so projections of future ice-sheet behaviour have to rely on bed topographies interpolated in a variety of ways between the direct measurements (Vaughan et al., 2006; Fretwell et al., 2013; Rignot et al., 2014; Millan et al., 2017; Morlighem et al., 2020). Over Thwaites Glacier, these interpolations have typically infilled areas of 15 km by 15 km between aerogeophysical flight lines, but 15 km between observations is much coarser than the resolution which Durand et al. (2011) suggest is desirable.

Where ice-penetrating radar surveys have been undertaken with sub-ice-thickness line spacing (e.g. Rutford Ice Stream – King et al., 2016; Pine Island Glacier – Bingham et al., 2017; Thwaites Glacier – Holschuh et al., 2020), they clearly identify details which are important for studying future ice-sheet behaviour that are not present in the interpolated bed topography products. In particular, imaged signatures in the bed often show some similarity to the much subtler topography of the ice surface above them. Theoretical studies based on linear perturbation theory (Gudmundsson, 2003; Raymond, 2005; Gudmundsson, 2008; Gudmundsson and Raymond, 2008; Raymond and Gudmundsson, 2009) have explored the relationship between the bed and the surface. The resulting relations can be used to infer bed properties from those of the surface but have only been applied twice to realistic settings: on 2D surface data from MacAyeal Ice Stream (Thorsteinsson et al., 2003) and on a 1D flow line from Rutford Ice Stream (Pralong and Gudmundsson, 2011). Both studies were undertaken in an era when surface elevation observations over Antarctica were of much lower quality and resolution than they are today.

Bed conditions such as geology, hydrology and sediment distribution also play a role in controlling ice flow and behaviour (Durand et al., 2011; Koellner et al., 2019) and are often poorly constrained. In many ice-sheet models, these bed conditions are combined into one parameter known as slipperiness, which is a measure of how easily the ice can slide over the topography (Rignot et al., 2011). Some seismic lines have been collected on Thwaites Glacier (Muto et al., 2019a, b), allowing a brief glimpse into the sediment distribution. Over the whole glacier, however, there are very few direct measurements of bed conditions which can be combined into slipperiness.

In this paper we exploit the relatively new availability of high-resolution surface elevation (REMA, ∼ 8 m; Howat et al., 2019) and velocity (NASA ITS_LIVE, ∼ 120 m resolution; Gardner et al., 2018) datasets. We apply linear perturbation theories to explore bed topography and slipperiness across the Thwaites Glacier catchment. We use a steady-state version of the shallow-ice-stream equations presented by Gudmundsson (2008) and compare the topography output from the inversion to radar grids and flight lines to assess its performance.

2 Methodology

2.1 Derivation of the steady-state shallow-ice-stream transfer functions

Gudmundsson (2008) derived a set of transfer functions which describe the relationship between the time-variant Fourier transforms of bed topography ($\widehat{b}$), bed slipperiness ($\widehat{c}$), surface topography ($\widehat{s}$) and horizontal components of surface velocity ($\widehat{u},\widehat{v}$). For the purposes considered here, this derivation can be simplified by considering the steady state from the beginning, removing the need to do a Laplace transform. Other than this, the derivation largely follows that of Gudmundsson (2008), but for clarity we state key assumptions and results here.

2.1.1 Response of flow to basal topography perturbations

Following Gudmundsson (2008) and working in a coordinate system tilted forward in the x direction by the mean surface slope, α, we start with the shallow-ice-stream equations of motion (MacAyeal, 1989),

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\partial }_x(\mathrm{4}h\mathit{\eta }{\partial }_{x}u& +\mathrm{2}h\mathit{\eta }{\partial }_{y}v)+{\partial }_y(h\mathit{\eta }({\partial }_{x}v+{\partial }_{y}u))\\ & -(u/c{)}^{\mathrm{1}/m}=\mathit{\rho }gh{\partial }_x(s\left)\text{cos}\right(\mathit{\alpha })-\mathit{\rho }gh\text{sin}(\mathit{\alpha }),\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\partial }_y(\mathrm{4}h\mathit{\eta }{\partial }_{y}v& +\mathrm{2}h\mathit{\eta }{\partial }_{x}u)+{\partial }_x(h\mathit{\eta }({\partial }_{y}u+{\partial }_{x}v))\\ & -(v/c{)}^{\mathrm{1}/m}=\mathit{\rho }gh{\partial }_y(s\left)\text{cos}\right(\mathit{\alpha }),\end{array}\end{array}$$

where u and v are the depth-independent velocity components in the x and y directions, respectively; h is the ice thickness; η is the effective ice viscosity; c is the basal slipperiness; m is a sliding-law parameter; ρ is the ice density; g is the acceleration due to gravity; s is the ice surface elevation; and α is the mean ice surface slope in the x direction.

Assuming that ice is a linear viscous medium (n=1) and that there is a non-linear sliding law (m>0), then the shallow-ice-stream equations can be linearised and solved analytically. We consider the spatial response to a small perturbation in basal topography, b, linearising around a reference model ($\stackrel{\mathrm{‾}}{h},\stackrel{\mathrm{‾}}{s},\stackrel{\mathrm{‾}}{b},\stackrel{\mathrm{‾}}{u},\stackrel{\mathrm{‾}}{v},\stackrel{\mathrm{‾}}{c}$) with $h=\stackrel{\mathrm{‾}}{h}+\mathrm{\Delta }h$, $s=\stackrel{\mathrm{‾}}{s}+\mathrm{\Delta }s$, $b=\stackrel{\mathrm{‾}}{b}+\mathrm{\Delta }b$, $u=\stackrel{\mathrm{‾}}{u}+\mathrm{\Delta }u$, v=Δv, w=Δw and $c=\stackrel{\mathrm{‾}}{c}$. The zero-order solutions are spatially constant, representing uniform flow down an inclined plane.

We, however, are interested in the first-order momentum balance equations:

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{xx}^{\mathrm{2}}\mathrm{\Delta }u+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{xy}^{\mathrm{2}}\mathrm{\Delta }v& +\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{yy}^{\mathrm{2}}\mathrm{\Delta }u-\mathit{\gamma }\mathrm{\Delta }u\\ & =\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\mathit{\alpha }{\partial }_x\mathrm{\Delta }s-\mathit{\rho }g\text{sin}\mathit{\alpha }\mathrm{\Delta }h,\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{yy}^{\mathrm{2}}\mathrm{\Delta }v+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{xy}^{\mathrm{2}}\mathrm{\Delta }u& +\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{xx}^{\mathrm{2}}\mathrm{\Delta }v-\mathit{\gamma }\mathrm{\Delta }v\\ & =\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\mathit{\alpha }{\partial }_y\mathrm{\Delta }s.\end{array}\end{array}$$

Also to the first order and importantly in the steady state, we have the following upper and lower kinematic boundary conditions:

$$\begin{array}{}\text{(5)}& {\displaystyle }\stackrel{\mathrm{‾}}{u}{\partial }_x\mathrm{\Delta }s-\mathrm{\Delta }w\left(s\right)=\mathrm{0},\text{(6)}& {\displaystyle }\stackrel{\mathrm{‾}}{u}{\partial }_x\mathrm{\Delta }b-\mathrm{\Delta }w\left(b\right)=\mathrm{0}.\end{array}$$

Various points about the validity of the steady-state assumption for Thwaites Glacier are raised in the discussion (Sect. 4).

All variables are then Fourier transformed with respect to the spatial variables x and y. Fourier-transformed variables are denoted with a circumflex ( $\widehat{}$ ). In the forward Fourier transform the wavenumbers in the x and y directions are denoted by k and l, respectively. This Fourier transform gives

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}\widehat{u}& +\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\widehat{v}+\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}\widehat{u}+\mathit{\gamma }\widehat{u}\\ & =ik\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\left(\mathit{\alpha }\right)\widehat{s}+\mathit{\rho }g\text{sin}\left(\mathit{\alpha }\right)\widehat{h},\end{array}\text{(8)}& {\displaystyle }\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}\widehat{v}+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\widehat{u}+\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}\widehat{v}+\mathit{\gamma }\widehat{v}=\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\left(\mathit{\alpha }\right)il\widehat{s},\text{(9)}& {\displaystyle }\widehat{w}\left(\stackrel{\mathrm{‾}}{s}\right)=-i\stackrel{\mathrm{‾}}{u}k\widehat{s},\text{(10)}& {\displaystyle }\widehat{w}\left(\stackrel{\mathrm{‾}}{b}\right)=-i\stackrel{\mathrm{‾}}{u}k\widehat{b},\end{array}$$

where $\widehat{h}=\widehat{s}-\widehat{b}$.

From depth integration of the Fourier-transformed incompressibility condition $w_z+u_x+v_y=\mathrm{0}$, we have

$$\begin{array}{}\text{(11)}& i\stackrel{\mathrm{‾}}{h}(k\widehat{u}+l\widehat{v})=\widehat{w}\left(\stackrel{\mathrm{‾}}{s}\right)-\widehat{w}\left(\stackrel{\mathrm{‾}}{b}\right),\end{array}$$

which, along with the steady-state boundary conditions, yields

$$\begin{array}{}\text{(12)}& i\stackrel{\mathrm{‾}}{h}(k\widehat{u}+l\widehat{v})=-ik\stackrel{\mathrm{‾}}{u}\widehat{s}+ik\stackrel{\mathrm{‾}}{u}\widehat{b}.\end{array}$$

Equations (7), (8) and (12) form a linear system of equations in $\widehat{s}$, $\widehat{u}$, $\widehat{v}$ and $\widehat{b}$ which can be solved algebraically (see Appendix A), leading to the steady-state transfer functions

$$\begin{array}{}\text{(13)}& {\displaystyle }{T}_{\text{sb}}(k,l)& {\displaystyle }={\displaystyle \frac{\widehat{s}}{\widehat{b}}}={\displaystyle \frac{\mathrm{i}k(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}})}{p\mathit{\xi }}},\text{(14)}& {\displaystyle }{T}_{\text{ub}}(k,l)& {\displaystyle }={\displaystyle \frac{\widehat{u}}{\widehat{b}}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }(l^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}-k^{\mathrm{2}}\stackrel{\mathrm{‾}}{u})}{\mathit{\xi }\mathit{\nu }p}},\text{(15)}& {\displaystyle }{T}_{\text{vb}}(k,l)& {\displaystyle }={\displaystyle \frac{\widehat{v}}{\widehat{b}}}={\displaystyle \frac{kl{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }({\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }\stackrel{\mathrm{‾}}{u})}{\mathit{\xi }\mathit{\nu }p}},\end{array}$$

which represent the ratio of the Fourier components of the surface to the Fourier components of the bed as a function of wavenumber. The following abbreviations are used for simplicity in the derivation: $\mathit{\xi }=\mathit{\gamma }+\mathrm{4}\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}\mathit{\eta }$, $\mathit{\gamma }={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}^{\mathrm{1}-m}}{m\stackrel{\mathrm{‾}}{c}}}$, $j^{\mathrm{2}}=k^{\mathrm{2}}+l^{\mathrm{2}}$, ${\mathit{\tau }}_{\mathrm{d}}=\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{sin}\mathit{\alpha }$, $p={\displaystyle \frac{\mathrm{i}}{t_{\mathrm{p}}}}-{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{r}}}}$, ${\displaystyle \frac{\mathrm{1}}{t_{\mathrm{p}}}}=k(\stackrel{\mathrm{‾}}{u}+{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\mathit{\xi }}})$, ${\displaystyle \frac{\mathrm{1}}{t_{\mathrm{r}}}}={\displaystyle \frac{j^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{h}\text{cot}\mathit{\alpha }}{\mathit{\xi }}}$ and $\mathit{\nu }=\mathit{\gamma }+\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}\mathit{\eta }$.

2.1.2 Response of flow to basal slipperiness perturbation

Starting once again with the shallow-ice-stream equations (Eqs. 1 and 2; MacAyeal, 1989), this time we consider the response to a small perturbation in basal slipperiness, c, linearising with $h=\stackrel{\mathrm{‾}}{h}+\mathrm{\Delta }h=\stackrel{\mathrm{‾}}{h}+\mathrm{\Delta }s$, $s=\stackrel{\mathrm{‾}}{s}+\mathrm{\Delta }s$, $b=\stackrel{\mathrm{‾}}{b}$, $u=\stackrel{\mathrm{‾}}{u}+\mathrm{\Delta }u$, v=Δv, w=Δw and $c=\stackrel{\mathrm{‾}}{c}(\mathrm{1}+\mathrm{\Delta }c)$ where Δc is the fractional slipperiness.

This gives the following first-order momentum balance equations:

$$\begin{array}{}\text{(16)}& {\displaystyle }\begin{array}{rl}\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{xx}^{\mathrm{2}}\mathrm{\Delta }u& +\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{xy}^{\mathrm{2}}\mathrm{\Delta }v+\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{yy}^{\mathrm{2}}\mathrm{\Delta }u-\mathit{\gamma }\mathrm{\Delta }u\\ & =\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\mathit{\alpha }{\partial }_x\mathrm{\Delta }s-\mathit{\rho }g\text{sin}\mathit{\alpha }\mathrm{\Delta }s-\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\mathrm{\Delta }c,\end{array}\text{(17)}& {\displaystyle }\begin{array}{rl}\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{yy}^{\mathrm{2}}\mathrm{\Delta }v& +\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{xy}^{\mathrm{2}}\mathrm{\Delta }u+\mathit{\eta }\stackrel{\mathrm{‾}}{h}{\partial }_{xx}^{\mathrm{2}}\mathrm{\Delta }v-\mathit{\gamma }\mathrm{\Delta }v\\ & =\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\mathit{\alpha }{\partial }_y\mathrm{\Delta }s.\end{array}\end{array}$$

Fourier transforming with respect to the spatial variables x and y gives

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{rl}\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}\widehat{u}& +\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\widehat{v}+\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}\widehat{u}+\mathit{\gamma }\widehat{u}\\ & =\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\mathit{\alpha }ik\widehat{s}+\mathit{\rho }g\text{sin}\mathit{\alpha }\widehat{s}+\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c},\end{array}\text{(19)}& {\displaystyle }\begin{array}{rl}\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}\widehat{v}& +\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\widehat{u}+\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}\widehat{v}+\mathit{\gamma }\widehat{v}\\ & =\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\mathit{\alpha }il\widehat{s}.\end{array}\end{array}$$

As there is no bed topography perturbation, the steady-state boundary conditions become

$$\begin{array}{}\text{(20)}& i\stackrel{\mathrm{‾}}{h}(k\widehat{u}+l\widehat{v})=-ik\stackrel{\mathrm{‾}}{u}\widehat{s}.\end{array}$$

Equations (18)–(20) form a linear system of equations which can be solved using standard algebraic techniques (see Appendix B), leading to the steady-state transfer functions

$$\begin{array}{}\text{(21)}& {\displaystyle }{T}_{\text{sc}}(k,l)& {\displaystyle }={\displaystyle \frac{\widehat{s}}{\widehat{c}}}=-{\displaystyle \frac{\mathrm{i}k\stackrel{\mathrm{‾}}{h}\stackrel{\mathrm{‾}}{u}\mathit{\gamma }}{p\mathit{\xi }}},\text{(22)}& {\displaystyle }{T}_{\text{uc}}(k,l)& {\displaystyle }={\displaystyle \frac{\widehat{u}}{\widehat{c}}}={\displaystyle \frac{\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\left(\right(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu }\left)\right(ik\stackrel{\mathrm{‾}}{u})-l^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\stackrel{\mathrm{‾}}{h})}{\mathit{\xi }\mathit{\nu }p}},\text{(23)}& {\displaystyle }{T}_{\text{vc}}(k,l)& {\displaystyle }={\displaystyle \frac{\widehat{v}}{\widehat{c}}}={\displaystyle \frac{kl\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\stackrel{\mathrm{‾}}{h}({\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }-\mathrm{3}i\mathit{\eta }\stackrel{\mathrm{‾}}{u}k)}{\mathit{\xi }\mathit{\nu }p}},\end{array}$$

which represent the ratio of variability in the Fourier components of the surface to variability in the Fourier components of the slipperiness.

Note that the transfer functions T_uc and T_vc are not the same as the steady-state versions of the transfer functions published in Gudmundsson (2008) as there is a typographic error in their paper. However, when plotted graphically, they can be used to reproduce the figures in that paper.

2.1.3 Non-dimensionalisation

These transfer functions can also be considered in a non-dimensional form, allowing us to make more general statements about the behaviour of the system in terms of key variables, such as ice thickness as the characteristic length scale. For this purpose the same scalings as used in Gudmundsson (2003, 2008) are employed. All spatial scales are in units of mean ice thickness ($\stackrel{\mathrm{‾}}{h}$), and stress components are in units of driving stress (τ_d). Non-dimensional velocity components are in units of mean deformational velocity (${\stackrel{\mathrm{‾}}{u}}_{\mathrm{d}}$) where

$$\begin{array}{}\text{(24)}& {\stackrel{\mathrm{‾}}{u}}_{\mathrm{d}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{h}{\mathit{\tau }}_{\mathrm{d}}}{\mathrm{2}\mathit{\eta }}}.\end{array}$$

The scale for slipperiness is given by $\stackrel{\mathrm{‾}}{c}/\stackrel{\mathrm{‾}}{C}$, where $\stackrel{\mathrm{‾}}{c}$ is the mean dimensional slipperiness and $\stackrel{\mathrm{‾}}{C}$ is the mean non-dimensional slipperiness. $\stackrel{\mathrm{‾}}{C}$ depends on not only slipperiness but also viscosity and thickness (through its dependence on u_s and τ_d). From Gudmundsson (2008) we know that $\stackrel{\mathrm{‾}}{c}/\stackrel{\mathrm{‾}}{C}={\stackrel{\mathrm{‾}}{u}}_{\mathrm{d}}/{\mathit{\tau }}_{\mathrm{d}}^m$ and we also have $\stackrel{\mathrm{‾}}{C}={\stackrel{\mathrm{‾}}{u}}_{\mathrm{b}}/{\stackrel{\mathrm{‾}}{u}}_{\mathrm{d}}$. The ice surface velocity is the sum of the deformational velocity and the basal velocity, such that ${\stackrel{\mathrm{‾}}{u}}_{\mathrm{s}}={\stackrel{\mathrm{‾}}{u}}_{\mathrm{d}}+{\stackrel{\mathrm{‾}}{u}}_{\mathrm{b}}$. With some simple algebra, we can therefore express the scale for slipperiness in terms of the surface velocity ${\stackrel{\mathrm{‾}}{u}}_{\mathrm{s}}$, which is a known quantity,

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\stackrel{\mathrm{‾}}{c}}{\stackrel{\mathrm{‾}}{C}}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{s}}}{{\mathit{\tau }}_{\mathrm{d}}^m(\stackrel{\mathrm{‾}}{C}+\mathrm{1})}}.\end{array}$$

The non-dimensional form of the equations is obtained using the substitutions $\stackrel{\mathrm{‾}}{c}\mapsto \stackrel{\mathrm{‾}}{C}$, $\mathit{\eta }\mapsto \mathrm{1}/\mathrm{2}$, $\stackrel{\mathrm{‾}}{h}\mapsto \mathrm{1}$, $\stackrel{\mathrm{‾}}{u}\mapsto \stackrel{\mathrm{‾}}{C}$, $\mathit{\gamma }\mapsto (m\stackrel{\mathrm{‾}}{C}{)}^{-\mathrm{1}}$ and ${\mathit{\tau }}_{\mathrm{d}}=\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{sin}\mathit{\alpha }\mapsto \mathrm{1}$. The non-dimensional transfer functions can be found in the Supplement and are also shown graphically there. Non-dimensionalised parameters are represented by capital letters (B, C, S, U, V, T_SB, T_UB, T_VB, T_SC, T_UC and T_UV).

2.2 The inverse problem

The non-dimensional transfer functions (T_SB, T_UB, T_VB, T_SC, T_UC and T_VC) describe the relationship between the Fourier transforms of non-dimensionalised bed topography ($\widehat{B}$), bed slipperiness ($\widehat{C}$), surface topography ($\widehat{S}$) and surface velocity ($\widehat{U},\widehat{V}$), as functions of the wavenumbers k and l. If the bed topography and slipperiness are known, then surface topography and velocity components are given by the forward model:

$$\begin{array}{}\text{(26)}& {\displaystyle }\widehat{S}(k,l)=T_{\text{SB}}(k,l)\widehat{B}(k,l)+T_{\text{SC}}(k,l)\widehat{C}(k,l),\text{(27)}& {\displaystyle }\widehat{U}(k,l)=T_{\text{UB}}(k,l)\widehat{B}(k,l)+T_{\text{UC}}(k,l)\widehat{C}(k,l),\text{(28)}& {\displaystyle }\widehat{V}(k,l)=T_{\text{VB}}(k,l)\widehat{B}(k,l)+T_{\text{VC}}(k,l)\widehat{C}(k,l).\end{array}$$

For each set of wavenumbers in Fourier space (k and l), we have three known variables ($\widehat{S}(k,l)$, $\widehat{U}(k,l)$ and $\widehat{V}(k,l)$) and two unknowns ($\widehat{B}(k,l)$ and $\widehat{C}(k,l)$), so the system is over-determined. We can therefore solve these equations independently for each wavenumber component of non-dimensional bed topography and slipperiness using a weighted least-squares inversion of Eqs. (26)–(28). Short-wavelength features and features aligned with ice flow are problematic because they cause flow disturbances in the ice which do not reach the surface in a measurable way, and so they can not be inverted from the surface data. We include a filtering parameter (p_filt) to remove these problematic wavelength components such that increasing the filtering parameter (p_filt) removes progressively longer wavelength features. This method was first used by Thorsteinsson et al. (2003) in their study of MacAyeal Ice Stream (formerly Ice Stream E). The equations which solve Eqs. (26)–(28) are therefore not repeated here in the main text but are given in notation consistent with this paper in Appendix C.

2.3 Synthetic tests

Synthetic tests allow us to explore which bed features can and can not be resolved using this inversion method. First we create a synthetic bed topography (b) and slipperiness (c) on a 50 km by 50 km grid, with 120 m data spacing, purposely matching the data spacing of the ITS_LIVE velocity product (Gardner et al., 2018). We subtract the mean bed elevation, slope and slipperiness so that the bed varies about 0, and we taper the outer 5 km of the grid linearly to 0 at the edges. This reduces edge effects in the inversion because the Fourier transform requires a periodic domain. We tested a few other sensible tapering functions, including semi-sinusoids but observed negligible differences in the inversion output when compared to the linear function. We then non-dimensionalise the tapered bed using the length scales given in Sect. 2.1.3 and Fourier transform to obtain $\widehat{B}$ and $\widehat{C}$. The non-dimensional surface elevation ($\widehat{S}$) and the velocity components ($\widehat{U}$, $\widehat{V}$) are calculated using the forward model and dimensionalised using the length scales given in Sect. 2.1.3. To simulate measurement errors in the real surface data, we add random noise to the generated surface (s, u, v). This noise is white noise with a Gaussian low-pass filter applied in Fourier space to give it a non-random frequency distribution. We then taper, non-dimensionalise and Fourier transform the noise-added surface data. Finally, we invert for non-dimensional bed topography ($\widehat{B}$) and slipperiness ($\widehat{C}$) using the inversion procedure described in Sect. 2.2 and the Supplement. After dimensionalisation, the inverted bed can be compared to the synthetic bed to study the behaviour of the inversion.

2.3.1 Parameter value choices

When running synthetic tests, several model parameters can be varied, in addition to the synthetic bed topography and slipperiness. Following Gudmundsson (2008) and Thorsteinsson et al. (2003), the sliding-law constant was set to m=1 and the filtering parameter (Eq. C4) to $p_{\text{filt}}=-\mathrm{2}$. The mean ice thickness $\stackrel{\mathrm{‾}}{h}$, mean surface slope α and mean ice surface velocity $\stackrel{\mathrm{‾}}{u}$ depend on the region studied, but for these synthetic tests, they are set at $\stackrel{\mathrm{‾}}{h}$ = 2000 m, α = 0.002 and $\stackrel{\mathrm{‾}}{u}$ = 100 m s⁻¹, values thought to be appropriate for the Thwaites Glacier region (Gudmundsson, 2008; Howat et al., 2019; Morlighem et al., 2020; Gardner et al., 2018). After studying the results of the Thwaites Glacier inversion for a variety of values of $\stackrel{\mathrm{‾}}{C}$, the non-dimensional mean slipperiness, we set $\stackrel{\mathrm{‾}}{C}$ = 100. When applied to Eq. (25), these values give a dimensional mean slipperiness $\stackrel{\mathrm{‾}}{c}$ = 2.7 × 10⁻³ $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{1}}$. We set the weighting factors (Eq. C2) to be Σ_s=0.001, Σ_u=1 and Σ_v=1. This accounts for the mismatch in the relative magnitude of the non-dimensionalised surface elevation and velocity and means the least-squares inversion solves for all three factors equally.

Figure 1 The effect of orientation to flow direction on how well landforms (top row; created synthetically) can be resolved by the inversion. These tests are presented on a 50 km by 50 km grid, where in the inversion results (bottom row) the outer 5 km is greyed out to hide edge effects that will be neglected. In these simulations mean ice thickness $\stackrel{\mathrm{‾}}{h}$ = 2000 m, mean slipperiness $\stackrel{\mathrm{‾}}{C}$ is 100, surface slope α is 0.02, amplitude is 200 m and wavelength λ is 20 km. Landforms at an angle of less than 15^∘ to the ice flow direction are not well resolved because they fall in the null space of the inversion.

Figure 2 The effect of wavelength on how well landforms (top row; created synthetically) can be resolved by the inversion. These tests are presented on a 50 km by 50 km grid, where in the inversion results (bottom row) the outer 5 km is greyed out to hide edge effects that will be neglected. In these simulations mean ice thickness $\stackrel{\mathrm{‾}}{h}$ = 2000 m, mean slipperiness $\stackrel{\mathrm{‾}}{C}$ is 100, surface slope α is 0.02, amplitude is 200 m and angle to flow θ is 90^∘. Landforms with a wavelength of less than 2 km can not be well resolved due to the shallow-ice-stream approximation.

2.3.2 Resolution of bed forms

A two-dimensional Fourier transform decomposes an image into a weighted sum of two-dimensional sinusoidal basis functions. For this reason, all of our synthetic tests used sinusoidal bed topographies and slipperiness as these are the most illustrative of the capabilities of the inversion. Sinusoidal basis functions vary depending on three parameters: the wavenumbers in the x and y directions (k,l) and the weighting or amplitude of the sinusoid. However, rather than considering wavenumbers, it is more intuitive to consider the horizontal wavelength, λ, and angle, θ, to the direction of the flow, where $j^{\mathrm{2}}=l^{\mathrm{2}}+k^{\mathrm{2}}$, $\mathit{\lambda }/\stackrel{\mathrm{‾}}{h}=\mathrm{2}\mathit{\pi }/j$, k=jcos(θ) and l=jsin(θ).

Figure 3 The effect of amplitude on how well landforms and slipperiness (rows 1 and 3, respectively; created synthetically) can be resolved by the inversion. These tests are presented on a 50 km by 50 km grid, where in the inversion results (rows 2 and 4, respectively) the outer 5 km is greyed out to hide edge effects that will be neglected. In these simulations mean ice thickness $\stackrel{\mathrm{‾}}{h}$ is 2000 m, mean slipperiness $\stackrel{\mathrm{‾}}{C}$ is 100, surface slope α is 0.02, wavelength λ is 20 km and angle to flow θ is 60^∘. The limiting factor on how well small amplitudes can be resolved is errors in the ice surface data, which are simulated by adding noise to the synthetic surface generated from the synthetic bed. Noise is added with an amplitude of 2 m for the ice surface elevation (Howat et al., 2019), and an amplitude of 15 m s⁻¹ for the ice surface velocity (Gardner et al., 2018). Bedforms with an amplitude of less than 10 m, and slipperiness of less than around 1 × 10⁻⁴ $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{1}}$ are not well resolved.

Figures 1–3 show how well bed topography and slipperiness can be resolved by the inversion when the angle to flow, θ; wavelength, λ; and amplitude are varied. The amplitude of the noise added to the synthetically generated surface is ± 2 m for surface elevation and ± 15 m s⁻¹ for velocity components as these are at the upper limit of the errors for REMA (Howat et al., 2019) and ITS_LIVE (Gardner et al., 2018) in the Thwaites Glacier region.

These synthetic tests show that in this simple least-squares inversion, the bed can be well resolved if the angle to the flow is greater than 15^∘, the wavelength is more than 2000 m, and the amplitude is greater than 10 m for topography or 1.34 × 10⁻⁴ $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{1}}$ for variability in slipperiness. It is worth noting, however, that the resolution of wavelengths varies depending on the ice thickness, which is the non-dimensional scale factor for lengths. This means that ice thickness is directly proportional to the wavelength at which variations in bed topography and slipperiness should be resolvable.

2.4 Applying the inversion to real data

We now turn our attention to the methodology used to apply the synthetically tested inversion to real data, using the Thwaites Glacier catchment as our example.

Our base data for surface elevation and velocity were the REMA digital elevation model with 8 m resolution (Howat et al., 2019) and output from the NASA MEaSUREs ITS_LIVE project with 120 m resolution (Gardner et al., 2018), respectively. Based on the latter, we therefore inverted for bed topography and slipperiness at 120 m resolution. An estimate of the ice thickness in each 50 km by 50 km region is obtained from a 50 km averaged version of BedMachine Antarctica ice thickness (Morlighem et al., 2020), and this is the only prior information about ice thickness used in the inversion.

Figure 4 (a) Multiple overlapping grids are used in the inversion when applied to real data to allow the variability between solutions to be studied. (b) The effect of changing the number of overlapping grids and the fraction of the overlapping central region used to calculate the standard deviation of bed topography from the inversion in a roughly 60 km by 60 km region.

In the synthetic tests discussed above, the inversion was run over a single 50 km by 50 km grid, with the outer 10 % of the grid discarded to reduce edge effects introduced during the Fourier transform. To look at the whole Thwaites catchment, we could use a set of adjacent 50 km by 50 km grids, but instead we chose to use more densely distributed grids which overlap. For each grid point, we calculate nine different ((but overlapping) inverted beds and then the mean bed topography and standard deviation. The standard deviation is not a measure of the error, but since the main approximation in the physics is the linearisation, we interpret the standard deviation to be a measure of non-linearity. In each of the overlapping grids, we use a set of zero-order parameters (such as average ice thickness), and because these zero-order parameters vary between the grids, the linearisation is different and the resulting beds are also different. However, inappropriate application of the shallow-ice-stream approximation or edge effects could also be influencing this. Figure 4a summarises this methodology. The bed topography and slipperiness results presented here are the grid-point by grid-point means of nine overlapping grids where each overlapping region is 1.67 km by 1.67 km. These values were chosen following tests on a small region of the real data (Fig. 4b).

When applied to real surface elevation and velocity data, this method generates four products: the mean and standard deviation of bed topography and the mean and standard deviation of bed slipperiness. The standard deviation is not a measure of the error in the bed topography or bed slipperiness and should not be interpreted as such.

Figure 5 Inversion outputs across a 160 km by 280 km region of Thwaites Glacier (location shown in panels f and g). The compass directions shown are accurate at the centre of the region (red dot in panel g) but may vary by up to 10^∘ across the region due to the polar stereographic projection. (a) Bed topography from BedMachine Antarctica (Morlighem et al., 2020) and (b) bed topography, (c) standard deviation of bed topography, (d) bed slipperiness and (e) standard deviation of bed slipperiness from our inversions of REMA (Howat et al., 2019) and ITS_LIVE (Gardner et al., 2018) at 120 m resolution. Standard deviations are a measure of variability between overlapped patches and should not be interpreted as a measure of the error in the bed topography or bed slipperiness. Black rectangles depict regions of pre-existing highly resolved bed topography (Holschuh et al., 2020) examined in Fig. 6. The four subglacial lakes observed from surface altimetry changes by Smith et al. (2017) are also outlined in black.

3 Results for Thwaites Glacier

Figure 5 shows the bed conditions we inverted from REMA (Howat et al., 2019) and ITS_LIVE (Gardner et al., 2018) over a 280 km by 160 km region of the main trunk of Thwaites Glacier.

The bed topography product from the inversion is shown in Fig. 5b. On the basin scale, the main basal topographic features identified by the inversion are several sets of parallel ridges which are oriented perpendicularly to the direction of ice flow and smooth basins in between these ridges. The location of these ridges matches well with the BedMachine Antarctica bed (Fig. 5a), particularly around the subglacial lakes, which appear to be between successive sets of ridges. The smoother topography in the basins between ridges is reflected in the inversion, particularly in the basin to the east of the Upper Thwaites region (Basin Y, Fig. 5). Many smaller hills also match BedMachine Antarctica, such as those at the upstream (south) end of the Upper Thwaites radar grid. However, the inversion also generates some notable features that are not present in the BedMachine Antarctica bed, such as the north-eastern extent of the central ridge next to the most upstream subglacial lake (Ridge Z, Fig. 5b).

The standard deviation of the bed topography (Fig. 5c) represents the range in model outputs from overlapping grid squares which use different regions of the ice surface. As might be expected, this standard deviation is lower in the central trunk of the glacier where the topographic gradients are smaller. The standard deviation is higher at the edges of the glacier trunk where the gradient of the topography changes and the shallow-ice-stream approximation breaks down. Standard deviation is also high in the north-west part of the inversion where some of the input surface is the surface of the floating ice shelf, which is subject to different physical processes than the grounded ice. The inversion does not produce results for the north-west corner as there are gaps in the input data where the surface sampled is not ice but open ocean.

The bed slipperiness product from the inversion is shown in Fig. 5d. The pattern in the slipperiness output is similar to the topography, with a dominant east-to-west lineation, although it is slightly difficult to make out due to the strong underlying slipperiness variation across the Thwaites Glacier region. This directional trend in slipperiness is also observed in the slipperiness output of other inversions in the Thwaites Glacier region (Barnes et al., 2021). There is less variability in slipperiness in the basins where the topography is smoother. In contrast the ridges have much more variable slipperiness. Like the standard deviation of bed topography, the standard deviation of bed slipperiness (Fig. 5e) is highest around the edges of the central trunk where there are higher topographic gradients.

Figure 6 (a) Bed topography acquired at sub-ice-thickness resolution by swath radar (Holschuh et al., 2020); (b) our inverted bed topography with mean non-dimensional slipperiness $\stackrel{\mathrm{‾}}{C}$ = 100; and (c) our inverted bed topography with mean slipperiness $\stackrel{\mathrm{‾}}{C}$ = 150, for the site labelled “Upper Thwaites” in Fig. 5a. Panels (d–f) show equivalent products for the site labelled “Lower Thwaites” in Fig. 5a. The results of a simple linear regression between the swath radar bed and the inverted bed are also given, with r being the regression coefficient and slope the gradient of the line of best fit. A slope of 1 means the amplitude of the inverted bed matches the amplitude of the swath radar bed.

Figure 7 (a) Comparative plot of inverted bed topography with mean non-dimensional slipperiness $\stackrel{\mathrm{‾}}{C}$ = 100 with along-flow radar-sounded bed topography. Bed topography is given in three forms: unfiltered bed picks from the 2019/20 airborne surveys (Jordan and Robinson, 2021); a version of the same filtered to 2 km wavelengths to be more representative of the detail we might expect to image in our inversion; the bed profile extracted from BedMachine Antarctica (Morlighem et al., 2020). The envelope around the inverted bed topography shows plus or minus 1 standard deviation. Standard deviations are a measure of variability between overlapped patches and should not be interpreted as a measure of the error in the bed topography or bed slipperiness. The correlation coefficients (r) and slopes given are the results of a linear regression between the inverted bed or the BedMachine Antarctica bed and the filtered radar bed. (b) Profile location within the inverted grid. (c) As for panel (a) but with mean non-dimensional slipperiness $\stackrel{\mathrm{‾}}{C}$ = 150.

Figure 8 (a) Comparative plot of inverted bed topography with mean non-dimensional slipperiness $\stackrel{\mathrm{‾}}{C}$ = 100 with across-flow radar-sounded bed topography. Key as for Fig. 7.

3.1 Comparison to radar data

In order to assess the success of the inversion, the bed topography output can be compared to existing radar data. Figure 6 shows a comparison between the inverted bed topography and bed topography sounded by swath radar at sub-ice-thickness resolution across two 20 km by 40 km regions (Holschuh et al., 2020). The inverted bed shows a good match to the swath-radar-imaged bed at larger scales, picking out the locations of all the main hills and valleys. There is a better match for the Upper Thwaites region than the Lower Thwaites region, and the fact that the inversion detects the channel between the ridges in the downstream (left) part of Upper Thwaites is particularly encouraging.

We further compare the inverted bed topography with the bed topography sounded by swath radar along ice flow (Fig. 7) and by airborne radar across ice flow (Fig. 8) over Thwaites Glacier in the 2019/20 field season (Jordan and Robinson, 2021). Further comparisons to more radar flight lines collected in the same field season can be seen in the Supplement. These figures also demonstrate that the inversion performs well in detecting the main hills and valleys and also highlight that the amplitude of topographic features are not always resolved correctly. This is likely due to variability in the local mean slipperiness away from the imposed global value of non-dimensional slipperiness $\stackrel{\mathrm{‾}}{C}=\mathrm{100}$. If a non-dimensional slipperiness $\stackrel{\mathrm{‾}}{C}=\mathrm{150}$ is imposed (Figs. 6, 7c and 8c), then the amplitudes of the inverted topography are reduced. Sometimes there is also an offset between the inverted and radar-sounded beds, caused by using a 50 km averaged version of the BedMachine ice thickness as the inversion ice thickness, rather than more detailed prior information.

4 Discussion

Our results demonstrate significant promise for being able to invert for bed topography across parts of Antarctica and other polar regions from surface elevation and velocity datasets. Comparisons with existing radar data available from Thwaites Glacier suggest that within the central trunk of the glacier, the bed features identified by the inversion are normally in the correct locations but are not always centred around the correct depth. These average depth differences are primarily due to the mean ice thickness used in the inversion, which is a 50 km averaged version of the BedMachine Antarctica ice thickness (Morlighem et al., 2020). For regions where there are radar flight lines and grids, these radar observations could be used instead of the averaged ice thickness to ensure that the bed depth is correct. For regions where there are very few existing radar data, this method has the potential in the future to identify obstacles to flow which are significant enough to affect surface ice dynamics, even if there is uncertainty about the local or regional ice thickness.

The regions where the inverted bed deviates significantly from the topography picked from radar surveys are of particular interest in assessing the potential of our inversion. Differences between inverted topography and radar lines are likely to be due to physical processes which are not encapsulated by the shallow-ice-stream approximation. One place in which the shallow-ice-stream approximation is known to break down is where the mean slope of the bedrock becomes too steep (Gudmundsson, 2003; Le Meur et al., 2004). The effect of this can be observed around the edges of the central trunk of the glacier, where the topographic slope is steep and the match between the inverted bed and the radar lines is poor (Fig. 5). Gudmundsson (2003) derived a full-system non-hydrostatic momentum balance version of the transfer function used in this work. The full-system approach does not rely on the shallow-ice-stream approximation and should therefore perform better. Future work using this method is likely to incorporate these more complex equations.

A further consideration in comparing the inverted bed to real data is the steady-state assumption made when deriving the transfer functions. Without repeat radar measurements for Thwaites Glacier we can not be sure of the stability of the bed. If the bed beneath Thwaites Glacier is changing rapidly, as observed at Rutford Ice Stream (Smith et al., 2007), then the surface may not represent the current bed but some long-term average. However, observations at Pine Island Glacier (Brisbourne et al., 2017; Davies et al., 2018) suggest that the bed is not changing rapidly there, and it is possible that neighbouring Thwaites Glacier might be behaving similarly. Additionally, the erosion observed at Rutford Ice Stream does not significantly change the shape of the topography on the wavelengths resolved by this inversion. If drumlins or mega-scale glaciation lineations (MSGLs) were forming, we would not be able to detect them with this method, as landforms aligned to flow fall in the null space of the inversion.

The steady-state assumption applies not only to the bed but also to the ice surface. Ice surface lowering due to glacier thinning would also affect the steady-state assumption, but since generally the ice surface lowers in a relatively uniform way, this would not have a significant effect on the first-order variations in the ice surface or the results of the model. More significant would be changes in the ice surface due to the filling and draining of subglacial lakes, but these changes are normally fairly localised and would not propagate to the higher-wavelength Fourier components. For Thwaites Glacier, the location of subglacial lakes is relatively well known (Smith et al., 2017; Hoffman et al., 2020; Malczyk et al., 2020), and we predict troughs in these locations as expected. The ice surface also becomes more unstable closer to the grounding line, with increased crevassing which would affect the surface profile. However, since results in the region immediately adjacent to the grounding line are compromised by the different physics of the ice shelf anyway, this is not a significant concern. With these caveats, we therefore consider the steady-state assumption to be suitable for the purposes of this inversion.

If the steady-state assumption is valid, then the age of the datasets used in the inversion is not important. However, input data from different years or decades could also affect the steady-state assumption. The main surface expressions of known bed features appear to be fairly similar between REMA (Howat et al., 2019) (2008–2018) and the earlier Bamber DEM (Bamber et al., 2009) (2003–2008), supporting the validity of the steady-state assumption for Thwaites Glacier. However, we also note that non-steady-state changes in the ice surface may be the reason for some of the features we observe (such as Ridge Z, Fig. 5) in the inversion output which are not seen in the airborne flight lines.

As in any modelling study, it is important to explore the behaviour of the inversion when the parameters chosen are varied. In this inversion there are just four fixed parameters which are not derived from the input data: the sliding-law exponent, m; the filtering parameter, p_filt; the weighting factor, Σ_s; and the mean slipperiness, $\stackrel{\mathrm{‾}}{C}$. The filtering parameter p_filt controls which frequencies are suppressed in the inversion to avoid introducing singularities. Higher values of p_filt (closer to 0) will filter out higher frequencies (lower wavelengths), and so a value of $p_{\text{filt}}=-\mathrm{2}$ is chosen to filter out noisy short wavelengths while maintaining realistic bed features. The weighting factor Σ_s controls the balance in the inversion between the surface elevation and surface velocity data, with smaller values of Σ_s weighting the inversion towards the surface data. Varying p_filt and Σ_s for the inversion of the real surface data confirms the choice of values from the synthetic tests ($p_{\text{filt}}=-\mathrm{2}$, Σ_s=0.001) as sensible values which return the best match with real bed data (Figs. S3 and S4 in the Supplement).

There is less certainty over what is the most suitable value for the non-dimensional mean slipperiness $\stackrel{\mathrm{‾}}{C}$. Although we have some measurements of the bed properties of Thwaites Glacier from seismic lines, gravity and magnetic inversions (Diehl, 2008; Jordan et al., 2010; Muto et al., 2019a), these are spatially limited and it is not currently clear how these properties combine into slipperiness at the bed (Kyrke-Smith et al., 2017). If $\stackrel{\mathrm{‾}}{C}$ is higher, then the amplitude of bed variability in the inversion output falls. Given the geological variability likely to be associated with multiple rifted tectonic blocks (Dunham et al., 2020) and the sediments deposited in those rift basins (Muto et al., 2019a, b), it is unlikely that the mean slipperiness, $\stackrel{\mathrm{‾}}{C}$, is the same across the whole region modelled here. Modelling studies which compare the results of different inversion procedures show that slipperiness may be quite variable across the Thwaites Glacier catchment (Barnes et al., 2021). In addition we note that the trend is quite different from features observed in the inverted topography, showing that the slipperiness map is not a result of linear trade-offs with topography in the inversion solution. The three-dimensional radar grids (Holschuh et al., 2020) are both located within regions with more topographic variability, likely unlifted rift blocks. This may explain why a lower value of $\stackrel{\mathrm{‾}}{C}=\mathrm{100}$ gives the best match in these regions, whereas a higher value of $\stackrel{\mathrm{‾}}{C}=\mathrm{150}$ (more slippery) gives a better match for the radar lines which cover both lithologies. It may also be that the 3D grids, which contain both along- and across-flow variability, are more representative of the bed than the 2D radar lines. For this reason, the catchment-scale bed topography presented in Fig. 5 is from the inversion with $\stackrel{\mathrm{‾}}{C}=\mathrm{100}$.

This uncertainty in $\stackrel{\mathrm{‾}}{C}$ means that some prior radar information is useful in order to calibrate the inversion method to give the best results. However, changing $\stackrel{\mathrm{‾}}{C}$ only alters the amplitude, so even if there is no prior information, the inversion will still identify hills and troughs. More detailed analysis of seismic and gravity data alongside the results of the inversion could also reveal trends that could be useful when applying this technique elsewhere.

The high-resolution swath radar grids (Holschuh et al., 2020) presented here have already been included in BedMachine Antarctica. Using the radar grids which are currently available, we can not therefore explore how well this inversion method performs compared to BedMachine Antarctica (Morlighem et al., 2020) over a dense radar grid. Both techniques use ice surface elevation and velocity datasets. BedMachine Antarctica uses these datasets and the principles of mass conservation (or streamline diffusion in slow-moving areas) to interpolate between detailed prior information on ice thickness from existing radar measurements. In contrast, the inversion method only requires an estimate of average ice thickness for each 50 km by 50 km grid and uses the linear perturbation theory described. This 50 km averaged ice thickness is subtracted from the ice surface to provide a reference bed to which the inversion adds perturbations. However, even if a single ice-thickness ($\stackrel{\mathrm{‾}}{h}$) value is used for the entire catchment, then the inversion method will still identify the location of hills and troughs in the bed topography, although the amplitudes and absolute depths of these features may be affected. Since the mass conservation method used in BedMachine assesses the ice flow through a series of flux gates ideally constrained by topography, it is much more reliant on good-quality, closely spaced, ice-thickness measurements from radar systems.

.

Figure 9 The effect of ice-thickness resolution on the results of the inversion in the Lower Thwaites region, as explored in Fig. 6. (a) Inversion bed topography when the low-resolution (50 km average) ice-thickness values shown in (c) are used. (b) Inversion bed topography when the full-resolution (120 m) ice-thickness values from BedMachine Antarctica (Morlighem et al., 2020) (d) are used.

To explore the role that the 50 km averaged ice thickness plays in the results of the inversion, we computed a 50 km gridded version of the BedMachine Antarctica ice thickness (Fig. 9) and then carried out a new inversion over the Lower Thwaites region (where there is existing swath radar). In this alternate ice-thickness input, each 50 km by 50 km region contains only one ice-thickness value, which is the average over that region. The results of this re-run (Fig. 9) with reduced prior ice-thickness information show the same short-wavelength features as the results using the full ice-thickness input, illustrating that the inversion method presented here is not unduly influenced by the ice thickness derived from BedMachine Antarctica.

In Figs. 7 and 8, it appears that the method estimates shorter-wavelength topography more accurately than longer wavelengths. We demonstrate this in Fig. 10, which shows results after wavelengths greater than 50 km have been removed from all profiles. It is clear that the inversion identifies peaks and troughs in the bed, although the amplitude of these features is not always correct. Fourier components with a wavelength above 50 km mainly represent the prior ice-thickness information supplied to the inversion as this is the large-scale zero-order topography to which the first-order perturbations from the inversion are added. The greater match between the results of the inversion and the PASIN data (Jordan and Robinson, 2021) after this bandpass filter therefore provides further evidence that the ice thickness derived from BedMachine Antarctica does not influence the key results from the inversion.

Figure 10 The results of the inversion (orange) compared to PASIN radar flight lines (grey) and BedMachine Antarctica (light blue) for an along-flow profile and an across-flow profile (locations shown in panels e and f, respectively). Panels (b) and (d) show the bed profiles. Panels (a) and (c) show the results with any Fourier components between 40 and 50 km in wavelength progressively damped with a half-cosine filter and any Fourier components over 50 km in wavelength removed.

A comparison of the two methods over an area where radar data have not yet been incorporated into BedMachine would allow an assessment of the reliability of the two techniques and identification of any artificial bed features introduced by each. Since the two radar grids presented (Holschuh et al., 2020) were included in the derivation of BedMachine Antarctica, no independent test is possible until more radar grids are collected.

5 Conclusions

We present the method and results of an inversion of ice surface elevation and velocity for bed topography and slipperiness in the Thwaites Glacier region. Our method builds on the method used by Thorsteinsson et al. (2003) in their study of MacAyeal Ice Stream but is based on a steady-state linear perturbation analysis of the shallow-ice-stream equations (MacAyeal, 1989; Gudmundsson, 2003). Synthetic tests show that this method can resolve variability in bed topography and slipperiness on wavelengths greater than 1 ice thickness and at amplitudes of more than 10 m for topography or 1 × 10⁻⁴ $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{1}}$ for slipperiness, as long as the variability is not aligned with the ice flow direction. Comparison of the results of the inversion with radar grids and flight lines suggests that the inversion correctly identifies most short-wavelength (< 50 km horizontal) features in the bed, with the correlation coefficient of a linear regression between the inverted bed and the radar bed as high as r=0.93 along some flight lines. This method works best in the central trunk of the glacier, where the gradient of the long-wavelength topography is low and relatively constant. Mismatches between the inverted bed topography and radar measurements are probably due to one of three factors: an incorrect ice thickness for that region, an unusually sticky or slippery bed, or physical processes not accounted for by the steady-state linearised shallow-ice-stream approximation. Future work, including incorporating more local prior ice-thickness data from radar measurements and the non-hydrostatic transfer functions from Gudmundsson (2003), may help to reduce these mismatches. Overall, the inversion provides an additional tool for studying landforms in the bed beneath glaciers which have a significant impact on ice flow. It will be particularly useful in ice streams where radar flight lines are sparse and standard interpolation techniques struggle, potentially reducing uncertainties in modelling the future behaviour of those regions and their contributions to global sea-level rise.

Appendix A: Derivation of transfer functions from a topography perturbation

Starting from the shallow-ice-stream equations, linearising around a small perturbation in bed topography and taking the Fourier transform of the first-order equations, we have

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}\widehat{u}& +\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\widehat{v}+\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}\widehat{u}+\mathit{\gamma }\widehat{u}\\ & =\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\mathit{\alpha }ik\widehat{s}+\mathit{\rho }g\text{sin}\mathit{\alpha }\widehat{h},\end{array}\text{(8)}& {\displaystyle }\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}\widehat{v}+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\widehat{u}+\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}\widehat{v}+\mathit{\gamma }\widehat{v}=\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\mathit{\alpha }il\widehat{s},\text{(12)}& {\displaystyle }i\stackrel{\mathrm{‾}}{h}(k\widehat{u}+l\widehat{v})=-ik\stackrel{\mathrm{‾}}{u}\widehat{s}+ik\stackrel{\mathrm{‾}}{u}\widehat{b}.\end{array}$$

Equations (7), (8) and (12) form a linear system of equations in $\widehat{s}$, $\widehat{u}$, $\widehat{v}$ and $\widehat{b}$ which can be solved algebraically using standard techniques:

$$\begin{array}{rl}& \left[\begin{array}{cc}(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}+\mathit{\nu })& \left(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\right)\\ \left(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\right)& (\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\end{array}\right]\left[\begin{array}{c}\widehat{u}\\ \widehat{v}\end{array}\right]\\ & =\left[\begin{array}{c}(ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+\frac{{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}})\widehat{s}-\frac{{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}\widehat{b}\\ il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\end{array}\right].\end{array}$$

The determinant of the left-hand side of these equations is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}+\mathit{\nu })(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })-\left(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\right)\left(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\right)\\ {\displaystyle }=& {\displaystyle }\mathrm{9}{\mathit{\eta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{h}}^{\mathrm{2}}{k}^{\mathrm{2}}{l}^{\mathrm{2}}+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}\mathit{\nu }+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}\mathit{\nu }+{\mathit{\nu }}^{\mathrm{2}}-\mathrm{9}{\mathit{\eta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{h}}^{\mathrm{2}}{k}^{\mathrm{2}}{l}^{\mathrm{2}}\\ {\displaystyle }=& {\displaystyle }\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}\mathit{\nu }+(\stackrel{\mathrm{‾}}{h}\mathit{\eta }{j}^{\mathrm{2}}+\mathit{\gamma })\mathit{\nu }\\ {\displaystyle }=& {\displaystyle }(\mathrm{4}\stackrel{\mathrm{‾}}{h}\mathit{\eta }{j}^{\mathrm{2}}+\mathit{\gamma })\mathit{\nu }\\ {\displaystyle }=& {\displaystyle }\mathit{\xi }\mathit{\nu },\end{array}$$

where the following abbreviations have been used to make the algebra clearer to follow: $\mathit{\nu }=\stackrel{\mathrm{‾}}{h}\mathit{\eta }{j}^{\mathrm{2}}+\mathit{\gamma }$, $j^{\mathrm{2}}=l^{\mathrm{2}}+k^{\mathrm{2}}$ and $\mathit{\xi }=\mathrm{4}\stackrel{\mathrm{‾}}{h}\mathit{\eta }{j}^{\mathrm{2}}+\mathit{\gamma }$.

$$\begin{array}{rl}\left[\begin{array}{c}\widehat{u}\\ \widehat{v}\end{array}\right]=& {\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}\left[\begin{array}{cc}(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })& (-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl)\\ (-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl)& (\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}+\mathit{\nu })\end{array}\right]\\ & \times \left[\begin{array}{c}(ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+\frac{{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}})\widehat{s}-\frac{{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}\widehat{b}\\ il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\end{array}\right]\end{array}$$

Therefore we have

$$\begin{array}{c}{\displaystyle }\begin{array}{rl}\widehat{u}=& {\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}\left(\right(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu }\left)\right((ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}})\widehat{s}-{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}}\widehat{b})\\ & +(-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl)\left(il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\right)),\end{array}\\ {\displaystyle }\begin{array}{rl}\widehat{v}=& {\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}\left(\right(-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\left)\right((ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}})\widehat{s}-{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}}\widehat{b})\\ & +(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}+\mathit{\nu })\left(il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\right)),\end{array}\end{array}$$

which simplifies to

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}\widehat{u}={\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }\stackrel{\mathrm{‾}}{h}}}(& \mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\widehat{s}-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\widehat{b}+ik\mathit{\nu }\stackrel{\mathrm{‾}}{h}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\\ & +{\mathit{\tau }}_{\mathrm{d}}\mathit{\nu }\widehat{s}-{\mathit{\tau }}_{\mathrm{d}}\mathit{\nu }\widehat{b}),\end{array}\text{(A2)}& {\displaystyle }\widehat{v}={\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }\stackrel{\mathrm{‾}}{h}}}(-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl{\mathit{\tau }}_{\mathrm{d}}\widehat{s}+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl{\mathit{\tau }}_{\mathrm{d}}\widehat{b}+i\stackrel{\mathrm{‾}}{h}l\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}).\end{array}$$

We then have

$$\begin{array}{ll}{\displaystyle }i\stackrel{\mathrm{‾}}{h}(k\widehat{u}+l\widehat{v})=& {\displaystyle }{\displaystyle \frac{i}{\mathit{\xi }\mathit{\nu }}}\left(k\right(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\widehat{s}-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\widehat{b}+ik\mathit{\nu }\stackrel{\mathrm{‾}}{h}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\\ {\displaystyle }& {\displaystyle }+{\mathit{\tau }}_{\mathrm{d}}\mathit{\nu }\widehat{s}-{\mathit{\tau }}_{\mathrm{d}}\mathit{\nu }\widehat{b})+l(-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl{\mathit{\tau }}_{\mathrm{d}}\widehat{s}\\ {\displaystyle }& {\displaystyle }+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl{\mathit{\tau }}_{\mathrm{d}}\widehat{b}+i\stackrel{\mathrm{‾}}{h}l\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\left)\right)\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{i}{\mathit{\xi }\mathit{\nu }}}(i{k}^{\mathrm{2}}\mathit{\nu }\stackrel{\mathrm{‾}}{h}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}+k{\mathit{\tau }}_{\mathrm{d}}\mathit{\nu }\widehat{s}-k{\mathit{\tau }}_{\mathrm{d}}\mathit{\nu }\widehat{b}\\ {\displaystyle }& {\displaystyle }+i\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s})\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathit{\xi }}}(-j^{\mathrm{2}}\stackrel{\mathrm{‾}}{h}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}+ik{\mathit{\tau }}_{\mathrm{d}}\widehat{s}-ik{\mathit{\tau }}_{\mathrm{d}}\widehat{b}).\end{array}$$

In the steady state, the kinematic boundary condition is

$$\begin{array}{ll}{\displaystyle }i\stackrel{\mathrm{‾}}{h}(k\widehat{u}+l\widehat{v})& {\displaystyle }=\widehat{w}\left(s\right)-\widehat{w}\left(b\right)\\ {\displaystyle }& {\displaystyle }=-i\stackrel{\mathrm{‾}}{u}k\widehat{s}+i\stackrel{\mathrm{‾}}{u}k\widehat{b}.\end{array}$$

Substituting the expression from above gives

$${\displaystyle \frac{\mathrm{1}}{\mathit{\xi }}}(-j^{\mathrm{2}}\stackrel{\mathrm{‾}}{h}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}+ik{\mathit{\tau }}_{\mathrm{d}}\widehat{s}-ik{\mathit{\tau }}_{\mathrm{d}}\widehat{b})=-i\stackrel{\mathrm{‾}}{u}k\widehat{s}+i\stackrel{\mathrm{‾}}{u}k\widehat{b},$$

which can be rearranged as

$$\begin{array}{ll}{\displaystyle }(i\stackrel{\mathrm{‾}}{u}k+{\displaystyle \frac{ik{\mathit{\tau }}_{\mathrm{d}}}{\mathit{\xi }}}-{\displaystyle \frac{j^{\mathrm{2}}\stackrel{\mathrm{‾}}{h}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }}{\mathit{\xi }}})\widehat{s}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\xi }}}(ik{\mathit{\tau }}_{\mathrm{d}}+ik\stackrel{\mathrm{‾}}{u}\mathit{\xi })\widehat{b}\\ {\displaystyle }\left(ik\right(\stackrel{\mathrm{‾}}{u}+{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\mathit{\xi }}})-{\displaystyle \frac{j^{\mathrm{2}}\stackrel{\mathrm{‾}}{h}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }}{\mathit{\xi }}})\widehat{s}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\xi }}}\left(ik\right({\mathit{\tau }}_{\mathrm{d}}+\stackrel{\mathrm{‾}}{u}\mathit{\xi }\left)\right)\widehat{b}\\ {\displaystyle }\mathit{\xi }p\widehat{s}& {\displaystyle }=ik({\mathit{\tau }}_{\mathrm{d}}+\stackrel{\mathrm{‾}}{u}\mathit{\xi })\widehat{b}.\end{array}$$

In agreement with Gudmundsson (2008), this leads to the steady-state transfer function

$$\begin{array}{}\text{(13)}& T_{\text{sb}}(k,l,t)={\displaystyle \frac{\widehat{s}}{\widehat{b}}}={\displaystyle \frac{\mathrm{i}k(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}})}{p\mathit{\xi }}}.\end{array}$$

Expanding the expression for $\widehat{u}$ gives

$$\begin{array}{ll}{\displaystyle }\widehat{u}=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }\stackrel{\mathrm{‾}}{h}}}(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\widehat{s}-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\widehat{b}\\ \text{(A1 revisited)}& {\displaystyle }& {\displaystyle }+ik\mathit{\nu }\stackrel{\mathrm{‾}}{h}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}+{\mathit{\tau }}_{\mathrm{d}}\mathit{\nu }\widehat{s}-{\mathit{\tau }}_{\mathrm{d}}\mathit{\nu }\widehat{b}){\displaystyle }=& {\displaystyle }{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\mathit{\xi }\mathit{\nu }\stackrel{\mathrm{‾}}{h}}}\left(\stackrel{\mathrm{‾}}{h}\mathit{\nu }\right(ik\text{cot}\mathit{\alpha })\widehat{s}+(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\widehat{s}-(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu }\left)\widehat{b}\right)\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\mathit{\xi }\mathit{\nu }\stackrel{\mathrm{‾}}{h}}}\left(\stackrel{\mathrm{‾}}{h}\mathit{\nu }\right(ik\text{cot}\mathit{\alpha }\left)\right({\displaystyle \frac{\mathrm{i}k(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}})}{p\mathit{\xi }}})\widehat{b}+(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\\ {\displaystyle }& {\displaystyle }\times \left({\displaystyle \frac{\mathrm{i}k(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}})}{p\mathit{\xi }}}\right)\widehat{b}-(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\widehat{b})\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\mathit{\xi }\mathit{\nu }\stackrel{\mathrm{‾}}{h}p\mathit{\xi }}}\left(\stackrel{\mathrm{‾}}{h}\mathit{\nu }\right(ik\text{cot}\mathit{\alpha }\left)\right(\mathrm{i}k(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}}))+(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\\ {\displaystyle }& {\displaystyle }\times \left(\mathrm{i}k\right(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}}\left)\right)-(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })p\mathit{\xi })\widehat{b}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\mathit{\xi }\mathit{\nu }\stackrel{\mathrm{‾}}{h}p\mathit{\xi }}}\left(\stackrel{\mathrm{‾}}{h}\mathit{\nu }\right(ik\text{cot}\mathit{\alpha }\left)\right(\mathrm{i}k(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}}))+(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\\ {\displaystyle }& {\displaystyle }\times \left(\mathrm{i}k\right(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}}\left)\right)-(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\\ {\displaystyle }& {\displaystyle }\times \left(ik\right(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}})-j^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{h}\text{cot}\mathit{\alpha }))\widehat{b}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\mathit{\xi }\mathit{\nu }\stackrel{\mathrm{‾}}{h}p\mathit{\xi }}}\left(\stackrel{\mathrm{‾}}{h}\mathit{\nu }\right(ik\text{cot}\mathit{\alpha }\left)\right(\mathrm{i}k(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}}))-(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\\ {\displaystyle }& {\displaystyle }\times (-j^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{h}\text{cot}\mathit{\alpha }))\widehat{b}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(-k^{\mathrm{2}}\mathit{\nu }\stackrel{\mathrm{‾}}{u}\mathit{\xi }-\mathit{\nu }{k}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }{j}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}})\widehat{b}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}-k^{\mathrm{2}}\stackrel{\mathrm{‾}}{u}\mathit{\xi })\widehat{b}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(\mathit{\xi }{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}-k^{\mathrm{2}}\stackrel{\mathrm{‾}}{u}\mathit{\xi })\widehat{b},\\ {\displaystyle }\widehat{u}=& {\displaystyle }{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }}{\mathit{\xi }\mathit{\nu }p}}(l^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}-k^{\mathrm{2}}\stackrel{\mathrm{‾}}{u})\widehat{b},\end{array}$$

remembering that $\mathit{\xi }=\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}+\mathit{\nu }$.

In agreement with Gudmundsson (2008), this leads to the following steady-state transfer function:

$$\begin{array}{}\text{(14)}& {\displaystyle }{T}_{\text{ub}}(k,l,t){\displaystyle }={\displaystyle \frac{\widehat{u}}{\widehat{b}}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }(l^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}-k^{\mathrm{2}}\stackrel{\mathrm{‾}}{u})}{\mathit{\xi }\mathit{\nu }p}}.\end{array}$$

Expanding the expression for $\widehat{v}$ gives

$$\begin{array}{ll}{\displaystyle }\widehat{v}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }\stackrel{\mathrm{‾}}{h}}}(-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl{\mathit{\tau }}_{\mathrm{d}}\widehat{s}+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl{\mathit{\tau }}_{\mathrm{d}}\widehat{b}+i\stackrel{\mathrm{‾}}{h}l\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s})\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}l}{\mathit{\xi }\mathit{\nu }}}(-\mathrm{3}\mathit{\eta }k({\displaystyle \frac{\mathrm{i}k(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}})}{p\mathit{\xi }}})\widehat{b}+\mathrm{3}\mathit{\eta }k\widehat{b}\\ {\displaystyle }& {\displaystyle }+i\mathit{\nu }\text{cot}\mathit{\alpha }\left({\displaystyle \frac{\mathrm{i}k(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}})}{p\mathit{\xi }}}\right)\widehat{b})\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}l}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(-\mathrm{3}\mathit{\eta }k(\mathrm{i}k(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}}))+\mathrm{3}\mathit{\eta }kp\mathit{\xi }\\ {\displaystyle }& {\displaystyle }+i\mathit{\nu }\text{cot}\mathit{\alpha }\left(\mathrm{i}k\right(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}}\left)\right))\widehat{b}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}l}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(-\mathrm{3}\mathit{\eta }k(\mathrm{i}k(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}}))\\ {\displaystyle }& {\displaystyle }+\mathrm{3}\mathit{\eta }k\left(ik\right(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}})-j^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{h}\text{cot}\mathit{\alpha })\\ {\displaystyle }& {\displaystyle }+i\mathit{\nu }\text{cot}\mathit{\alpha }\left(\mathrm{i}k\right(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}}\left)\right))\widehat{b}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}lk}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(+\mathrm{3}\mathit{\eta }(-j^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{h}\text{cot}\mathit{\alpha })-\mathit{\nu }\text{cot}\mathit{\alpha }(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}}\left)\right)\widehat{b}\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}lk\text{cot}\mathit{\alpha }}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(-(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}+\mathit{\nu }){\mathit{\tau }}_{\mathrm{d}}-\mathit{\nu }\stackrel{\mathrm{‾}}{u}\mathit{\xi })\widehat{b},\\ {\displaystyle }\widehat{v}& {\displaystyle }={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}lk\text{cot}\mathit{\alpha }}{\mathit{\xi }\mathit{\nu }p}}({\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }\stackrel{\mathrm{‾}}{u})\widehat{b},\end{array}$$

remembering that $\mathit{\xi }=\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}+\mathit{\nu }$.

In agreement with Gudmundsson (2008), this leads to the following steady-state transfer function:

$$\begin{array}{}\text{(15)}& T_{\text{vb}}(k,l,t)={\displaystyle \frac{\widehat{v}}{\widehat{b}}}={\displaystyle \frac{kl{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }({\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }\stackrel{\mathrm{‾}}{u})}{\mathit{\xi }\mathit{\nu }p}}.\end{array}$$

Appendix B: Derivation of transfer functions from a slipperiness perturbation

Starting once again with the shallow-ice-stream equations (MacAyeal, 1989), this time we consider the response to a small perturbation in basal slipperiness, linearising and taking the Fourier transform to give

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{rl}\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}\widehat{u}& +\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\widehat{v}+\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}\widehat{u}+\mathit{\gamma }\widehat{u}\\ & =\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\mathit{\alpha }ik\widehat{s}+\mathit{\rho }g\text{sin}\mathit{\alpha }\widehat{s}+\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c},\end{array}\text{(19)}& {\displaystyle }\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}\widehat{v}+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\widehat{u}+\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}\widehat{v}+\mathit{\gamma }\widehat{v}=\mathit{\rho }g\stackrel{\mathrm{‾}}{h}\text{cos}\mathit{\alpha }il\widehat{s},\text{(20)}& {\displaystyle }i\stackrel{\mathrm{‾}}{h}(k\widehat{u}+l\widehat{v})=-ik\stackrel{\mathrm{‾}}{u}\widehat{s}.\end{array}$$

Equations (18)–(20) form a linear system of equations which can be solved using standard algebraic techniques:

$$\begin{array}{rl}& \left[\begin{array}{cc}(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}+\mathit{\nu })& \left(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\right)\\ \left(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\right)& (\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\end{array}\right]\left[\begin{array}{c}\widehat{u}\\ \widehat{v}\end{array}\right]\\ & =\left[\begin{array}{c}(ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+\frac{{\mathit{\tau }}_{\mathrm{d}}}{h})\widehat{s}+\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c}\\ il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\end{array}\right].\end{array}$$

The determinant of the left-hand side is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}+\mathit{\nu })(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })-\left(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\right)\left(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\right)\\ {\displaystyle }=& {\displaystyle }\mathrm{9}{\mathit{\eta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{h}}^{\mathrm{2}}{k}^{\mathrm{2}}{l}^{\mathrm{2}}+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}\mathit{\nu }+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}\mathit{\nu }+{\mathit{\nu }}^{\mathrm{2}}-\mathrm{9}{\mathit{\eta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{h}}^{\mathrm{2}}{k}^{\mathrm{2}}{l}^{\mathrm{2}}\\ {\displaystyle }=& {\displaystyle }\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}\mathit{\nu }+(\stackrel{\mathrm{‾}}{h}\mathit{\eta }{j}^{\mathrm{2}}+\mathit{\gamma })\mathit{\nu }\\ {\displaystyle }=& {\displaystyle }(\mathrm{4}\stackrel{\mathrm{‾}}{h}\mathit{\eta }{j}^{\mathrm{2}}+\mathit{\gamma })\mathit{\nu }\\ {\displaystyle }=& {\displaystyle }\mathit{\xi }\mathit{\nu },\end{array}$$

where the following abbreviations have been used to make the algebra easier to follow: $j^{\mathrm{2}}=l^{\mathrm{2}}+k^{\mathrm{2}}$, $\mathit{\xi }=\mathrm{4}\stackrel{\mathrm{‾}}{h}\mathit{\eta }{j}^{\mathrm{2}}+\mathit{\gamma }$ and $\mathit{\nu }=\stackrel{\mathrm{‾}}{h}\mathit{\eta }{j}^{\mathrm{2}}+\mathit{\gamma }$.

$$\begin{array}{rl}\left[\begin{array}{c}\widehat{u}\\ \widehat{v}\end{array}\right]=& {\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}\left[\begin{array}{cc}(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })& (-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl)\\ (-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl)& (\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}+\mathit{\nu })\end{array}\right]\\ & \times \left[\begin{array}{c}(ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+\frac{{\mathit{\tau }}_{\mathrm{d}}}{h})\widehat{s}+\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c}\\ il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\end{array}\right]\end{array}$$

Therefore we have

$$\begin{array}{c}{\displaystyle }\begin{array}{rl}\widehat{u}=& {\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}\left(\right(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu }\left)\right((ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{h}})\widehat{s}+\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c})\\ & +(-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl)\left(il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\right)),\end{array}\\ {\displaystyle }\begin{array}{rl}\widehat{v}=& {\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}\left(\right(-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\left)\right((ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{h}})\widehat{s}+\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c})\\ & +(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{k}^{\mathrm{2}}+\mathit{\nu })\left(il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\right)),\end{array}\end{array}$$

which simplifies to

$$\begin{array}{}\text{(B1)}& {\displaystyle }\begin{array}{rl}\widehat{u}=& {\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}\left(\right(\mathrm{3}\mathit{\eta }{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+{\displaystyle \frac{\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}})\widehat{s}\\ & +(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c}),\end{array}\text{(B2)}& {\displaystyle }\widehat{v}={\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}(-\mathrm{3}\mathit{\eta }kl{\mathit{\tau }}_{\mathrm{d}}\widehat{s}+\mathit{\nu }il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c}).\end{array}$$

We then have

$$\begin{array}{ll}{\displaystyle }i\stackrel{\mathrm{‾}}{h}(k\widehat{u}+l\widehat{v})=& {\displaystyle }{\displaystyle \frac{i\stackrel{\mathrm{‾}}{h}}{\mathit{\xi }\mathit{\nu }}}\left(k\right((\mathrm{3}\mathit{\eta }{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+{\displaystyle \frac{\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}})\widehat{s}\\ {\displaystyle }& {\displaystyle }+(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c})\\ {\displaystyle }& {\displaystyle }+l(-\mathrm{3}\mathit{\eta }kl{\mathit{\tau }}_{\mathrm{d}}\widehat{s}+\mathit{\nu }il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c}))\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{i\stackrel{\mathrm{‾}}{h}}{\mathit{\xi }\mathit{\nu }}}\left(k\right((\mathit{\nu }ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+{\displaystyle \frac{\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}})\widehat{s}+\mathit{\nu }\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c})\\ {\displaystyle }& {\displaystyle }+l\left(\mathit{\nu }il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}\right))\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{i\stackrel{\mathrm{‾}}{h}}{\mathit{\xi }}}(i{k}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}+{\displaystyle \frac{k{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}}\widehat{s}+k\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c}\\ {\displaystyle }& {\displaystyle }+i{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s})\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{i}{\mathit{\xi }}}(i\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}+k{\mathit{\tau }}_{\mathrm{d}}\widehat{s}+k\stackrel{\mathrm{‾}}{h}\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c}).\end{array}$$

At a steady state, ${\widehat{s}}_{\mathrm{t}}=\mathrm{0}$, so we have the following boundary condition:

$$\begin{array}{ll}{\displaystyle }-ik\stackrel{\mathrm{‾}}{u}\widehat{s}& {\displaystyle }=i\stackrel{\mathrm{‾}}{h}(k\widehat{u}+l\widehat{v})\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{i}{\mathit{\xi }}}(ih{j}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}+k{\mathit{\tau }}_{\mathrm{d}}\widehat{s}+k\stackrel{\mathrm{‾}}{h}\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c}),\\ {\displaystyle }-ik\stackrel{\mathrm{‾}}{h}\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c}& {\displaystyle }=-h{j}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}+ik{\mathit{\tau }}_{\mathrm{d}}\widehat{s}+ik\stackrel{\mathrm{‾}}{u}\mathit{\xi }\widehat{s},\\ {\displaystyle \frac{\mathrm{1}}{\mathit{\xi }}}(-ik\stackrel{\mathrm{‾}}{h}\mathit{\gamma }\stackrel{\mathrm{‾}}{u})\widehat{c}& {\displaystyle }=(-{\displaystyle \frac{h{j}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }}{\mathit{\xi }}}+ik({\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{\mathit{\xi }}}+\stackrel{\mathrm{‾}}{u}\left)\right)\widehat{s},\\ {\displaystyle \frac{\mathrm{1}}{\mathit{\xi }}}(-ik\stackrel{\mathrm{‾}}{h}\mathit{\gamma }\stackrel{\mathrm{‾}}{u})\widehat{c}& {\displaystyle }=(-{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{r}}}}+{\displaystyle \frac{i}{t_{\mathrm{p}}}})\widehat{s}\\ {\displaystyle }(-ik\stackrel{\mathrm{‾}}{h}\mathit{\gamma }\stackrel{\mathrm{‾}}{u})\widehat{c}& {\displaystyle }=\mathit{\xi }p\widehat{s}.\end{array}$$

In agreement with Gudmundsson (2008) this leads to the following steady-state transfer function:

$$\begin{array}{}\text{(21)}& T_{\text{sc}}(k,l,t)={\displaystyle \frac{\widehat{s}(k,l,t)}{\widehat{c}(k,l)}}=-{\displaystyle \frac{\mathrm{i}k\stackrel{\mathrm{‾}}{h}\stackrel{\mathrm{‾}}{u}\mathit{\gamma }}{p\mathit{\xi }}}.\end{array}$$

Expanding the expression for $\widehat{u}$ (Eq. B1) gives

$$\begin{array}{ll}{\displaystyle }\widehat{u}=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}\left(\right(\mathrm{3}\mathit{\eta }{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+{\displaystyle \frac{\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}})\widehat{s}\\ {\displaystyle }& {\displaystyle }+(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c})\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}\left(\right(\mathrm{3}\mathit{\eta }{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }ik{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+{\displaystyle \frac{\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}}{\stackrel{\mathrm{‾}}{h}}}\left)\right(-{\displaystyle \frac{\mathrm{i}k\stackrel{\mathrm{‾}}{h}\stackrel{\mathrm{‾}}{u}\mathit{\gamma }}{p\mathit{\xi }}})\widehat{c}\\ {\displaystyle }& {\displaystyle }+(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c})\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathit{\gamma }\stackrel{\mathrm{‾}}{u}}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}\left(\right(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }i\stackrel{\mathrm{‾}}{h}k{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}\left)\right(-\mathrm{i}k)\\ {\displaystyle }& {\displaystyle }+(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })p\mathit{\xi })\widehat{c}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathit{\gamma }\stackrel{\mathrm{‾}}{u}}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}\left(\right(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }i\stackrel{\mathrm{‾}}{h}k{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }+\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}\left)\right(-\mathrm{i}k)\\ {\displaystyle }& {\displaystyle }+(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\left(ik\right(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}})-j^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{h}\text{cot}\mathit{\alpha }))\widehat{c}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathit{\gamma }\stackrel{\mathrm{‾}}{u}}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(\mathit{\nu }{k}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\stackrel{\mathrm{‾}}{h}+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}ik\stackrel{\mathrm{‾}}{u}\mathit{\xi }-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}{j}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{h}\text{cot}\mathit{\alpha }\\ {\displaystyle }& {\displaystyle }+\mathit{\nu }ik\stackrel{\mathrm{‾}}{u}\mathit{\xi }-\mathit{\nu }{j}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{h}\text{cot}\mathit{\alpha })\widehat{c}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathit{\gamma }\stackrel{\mathrm{‾}}{u}}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(-\mathit{\nu }{l}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\stackrel{\mathrm{‾}}{h}-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}{j}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{h}\text{cot}\mathit{\alpha }+\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}ik\stackrel{\mathrm{‾}}{u}\mathit{\xi }\\ {\displaystyle }& {\displaystyle }+\mathit{\nu }ik\stackrel{\mathrm{‾}}{u}\mathit{\xi })\widehat{c}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathit{\gamma }\stackrel{\mathrm{‾}}{u}}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(-(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}+\mathit{\nu })l^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\stackrel{\mathrm{‾}}{h}+(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu })\\ {\displaystyle }& {\displaystyle }\times \left(ik\stackrel{\mathrm{‾}}{u}\right)\mathit{\xi })\widehat{c},\\ {\displaystyle }\widehat{u}=& {\displaystyle }{\displaystyle \frac{\mathit{\gamma }\stackrel{\mathrm{‾}}{u}}{\mathit{\xi }\mathit{\nu }p}}(-l^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\stackrel{\mathrm{‾}}{h}+(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu }\left)\right(ik\stackrel{\mathrm{‾}}{u}\left)\right)\widehat{c},\end{array}$$

remembering that $\mathit{\xi }=\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}+\mathit{\nu }$, and leads to the following steady-state transfer function:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{T}_{\text{uc}}(k,l,t)& ={\displaystyle \frac{\widehat{u}}{\widehat{c}}}\\ & ={\displaystyle \frac{\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\left(\right(\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{l}^{\mathrm{2}}+\mathit{\nu }\left)\right(ik\stackrel{\mathrm{‾}}{u})-l^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\stackrel{\mathrm{‾}}{h})}{\mathit{\xi }\mathit{\nu }p}},\end{array}\end{array}$$

which is not as stated by Gudmundsson (2008).

Expanding the expression for $\widehat{v}$ (Eq. B2) gives

$$\begin{array}{ll}{\displaystyle }\widehat{v}=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}(-\mathrm{3}\mathit{\eta }kl{\mathit{\tau }}_{\mathrm{d}}\widehat{s}+\mathit{\nu }il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\widehat{s}-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c})\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathit{\xi }\mathit{\nu }}}(-\mathrm{3}\mathit{\eta }kl{\mathit{\tau }}_{\mathrm{d}}(-{\displaystyle \frac{\mathrm{i}k\stackrel{\mathrm{‾}}{h}\stackrel{\mathrm{‾}}{u}\mathit{\gamma }}{p\mathit{\xi }}})\widehat{c}+\mathit{\nu }il{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }(-{\displaystyle \frac{\mathrm{i}k\stackrel{\mathrm{‾}}{h}\stackrel{\mathrm{‾}}{u}\mathit{\gamma }}{p\mathit{\xi }}})\widehat{c}\\ {\displaystyle }& {\displaystyle }-\mathrm{3}\mathit{\eta }\stackrel{\mathrm{‾}}{h}kl\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\widehat{c})\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{kl\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\stackrel{\mathrm{‾}}{h}}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(\mathrm{3}i\mathit{\eta }k{\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }-\mathrm{3}\mathit{\eta }p\mathit{\xi })\widehat{c}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{kl\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\stackrel{\mathrm{‾}}{h}}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(\mathrm{3}i\mathit{\eta }k{\mathit{\tau }}_{\mathrm{d}}+\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }\\ {\displaystyle }& {\displaystyle }-\mathrm{3}\mathit{\eta }\left(ik\right(\stackrel{\mathrm{‾}}{u}\mathit{\xi }+{\mathit{\tau }}_{\mathrm{d}})-j^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{h}\text{cot}\mathit{\alpha }))\widehat{c}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{kl\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\stackrel{\mathrm{‾}}{h}}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}(\mathit{\nu }{\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }-\mathrm{3}\mathit{\eta }ik\stackrel{\mathrm{‾}}{u}\mathit{\xi }+\mathrm{3}\mathit{\eta }{j}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}\stackrel{\mathrm{‾}}{h}\text{cot}\mathit{\alpha })\widehat{c}\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{kl\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\stackrel{\mathrm{‾}}{h}}{\mathit{\xi }\mathit{\nu }p\mathit{\xi }}}\left(\right(\mathrm{3}\mathit{\eta }{j}^{\mathrm{2}}\stackrel{\mathrm{‾}}{h}+\mathit{\nu }){\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }-\mathrm{3}i\mathit{\eta }\stackrel{\mathrm{‾}}{u}k\mathit{\xi })\widehat{c},\\ {\displaystyle }\widehat{v}=& {\displaystyle }{\displaystyle \frac{kl\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\stackrel{\mathrm{‾}}{h}}{\mathit{\xi }\mathit{\nu }p}}({\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }-\mathrm{3}i\mathit{\eta }\stackrel{\mathrm{‾}}{u}k)\widehat{c},\end{array}$$

remembering that $\mathit{\xi }=\mathrm{4}\mathit{\eta }\stackrel{\mathrm{‾}}{h}{j}^{\mathrm{2}}+\mathit{\nu }$, which leads to the following steady-state transfer function:

$$\begin{array}{}\text{(23)}& T_{\text{vc}}(k,l,t)={\displaystyle \frac{\widehat{v}}{\widehat{c}}}={\displaystyle \frac{kl\mathit{\gamma }\stackrel{\mathrm{‾}}{u}\stackrel{\mathrm{‾}}{h}({\mathit{\tau }}_{\mathrm{d}}\text{cot}\mathit{\alpha }-\mathrm{3}i\mathit{\eta }\stackrel{\mathrm{‾}}{u}k)}{\mathit{\xi }\mathit{\nu }p}},\end{array}$$

which is also not as stated by Gudmundsson (2008).

Appendix C: The inverse problem

The transfer functions (T_sb, T_ub, T_vb, T_sc, T_uc and T_vc) describe the relationship between the Fourier transforms of bed topography ($\widehat{b}$), bed slipperiness ($\widehat{c}$), surface topography ($\widehat{s}$) and surface velocity ($\widehat{u},\widehat{v}$). If the bed topography and slipperiness are known, then surface topography and velocity components are given by the forward model:

$$\begin{array}{}\text{(26)}& {\displaystyle }\widehat{s}=T_{\text{sb}}\widehat{b}+T_{\text{sc}}\widehat{c},\text{(27)}& {\displaystyle }\widehat{u}=T_{\text{ub}}\widehat{b}+T_{\text{uc}}\widehat{c},\text{(28)}& {\displaystyle }\widehat{v}=T_{\text{vb}}\widehat{b}+T_{\text{vc}}\widehat{c}.\end{array}$$

Non-dimensionalised this gives

$$\begin{array}{}\text{(Eq. 26 non-dimensionalised)}& {\displaystyle }\widehat{S}=T_{\text{SB}}\widehat{B}+T_{\text{SC}}\widehat{C},\text{(Eq. 27 non-dimensionalised)}& {\displaystyle }\widehat{U}=T_{\text{UB}}\widehat{B}+T_{\text{UC}}\widehat{C},\text{(Eq. 28 non-dimensionalised)}& {\displaystyle }\widehat{V}=T_{\text{VB}}\widehat{B}+T_{\text{VC}}\widehat{C}.\end{array}$$

Since the system is over-determined, we can use a weighted least-squares inversion of Eqs. (26)–(28) to find the bed topography and slipperiness which are the most consistent with the ice surface. This is the same method used by Thorsteinsson et al. (2003) in their study of MacAyeal Ice Stream (formerly Ice Stream E) but is reproduced here in notation consistent with the rest of the equations we present.

In matrix form we have the following forward model:

$$\begin{array}{}\text{(C1)}& \mathbf{Y}=\mathbf{GX},\end{array}$$

where $\mathbf{Y}=\left[\begin{array}{c}\widehat{S}\\ \widehat{U}\\ \widehat{V}\end{array}\right]$, $\mathbf{G}=\left[\begin{array}{cc}{T}_{\text{SB}}& T_{\text{SC}}\\ T_{\text{UB}}& T_{\text{UC}}\\ T_{\text{VB}}& T_{\text{VC}}\end{array}\right]$, $\mathbf{X}=\left[\begin{array}{c}\widehat{B}\\ \widehat{C}\end{array}\right]$ and E = $\left[\begin{array}{ccc}{\sum }_S^{-\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\sum }_U^{-\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\sum }_V^{-\mathrm{2}}\end{array}\right]$.

A least-squares inversion gives

$$\begin{array}{ll}{\displaystyle }\mathbf{X}=& {\displaystyle }({\mathbf{G}}^H\mathbf{EG}{)}^{-\mathrm{1}}{\mathbf{G}}^H\mathbf{EY},\\ {\displaystyle }\left[\begin{array}{c}\widehat{B}\\ \widehat{C}\end{array}\right]=& {\displaystyle }\left(\left[\begin{array}{ccc}{T}_{\text{SB}}^{\ast }& T_{\text{UB}}^{\ast }& T_{\text{VB}}^{\ast }\\ T_{\text{SC}}^{\ast }& T_{\text{UC}}^{\ast }& T_{\text{VC}}^{\ast }\end{array}\right]\left[\begin{array}{ccc}{\sum }_S^{-\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\sum }_U^{-\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\sum }_V^{-\mathrm{2}}\end{array}\right]\right.\\ {\displaystyle }& {\displaystyle }{\left.\times \left[\begin{array}{cc}{T}_{\text{SB}}& T_{\text{SC}}\\ T_{\text{UB}}& T_{\text{UC}}\\ T_{\text{VB}}& T_{\text{VC}}\end{array}\right]\right)}^{-\mathrm{1}}\left[\begin{array}{ccc}{T}_{\text{SB}}^{\ast }& T_{\text{UB}}^{\ast }& T_{\text{VB}}^{\ast }\\ T_{\text{SC}}^{\ast }& T_{\text{UC}}^{\ast }& T_{\text{VC}}^{\ast }\end{array}\right]\\ \text{(C2)}& {\displaystyle }& {\displaystyle }\times \left[\begin{array}{ccc}{\sum }_S^{-\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\sum }_U^{-\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\sum }_V^{-\mathrm{2}}\end{array}\right]\left[\begin{array}{c}\widehat{S}\\ \widehat{U}\\ \widehat{V}\end{array}\right],\end{array}$$

where ^H is the Hermitian transpose and ^∗ is the complex conjugate.

For compactness we define the following:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }K=T_{\text{SB}}^{\ast }{T}_{\text{SC}}/\sum _S^{\mathrm{2}}+T_{\text{UB}}^{\ast }{T}_{\text{UC}}/\sum _U^{\mathrm{2}}+T_{\text{VB}}^{\ast }{T}_{\text{VC}}/\sum _V^{\mathrm{2}},\\ {\displaystyle }& {\displaystyle }L=T_{\text{SB}}^{\ast }{T}_{\text{SB}}/\sum _S^{\mathrm{2}}+T_{\text{UB}}^{\ast }{T}_{\text{UB}}/\sum _U^{\mathrm{2}}+T_{\text{VB}}^{\ast }{T}_{\text{VB}}/\sum _V^{\mathrm{2}},\\ {\displaystyle }& {\displaystyle }M=T_{\text{SC}}^{\ast }{T}_{\text{SC}}/\sum _S^{\mathrm{2}}+T_{\text{UC}}^{\ast }{T}_{\text{UC}}/\sum _U^{\mathrm{2}}+T_{\text{VC}}^{\ast }{T}_{\text{VC}}/\sum _V^{\mathrm{2}},\\ {\displaystyle }& {\displaystyle }{Y}_{\mathrm{b}}=\widehat{S}{T}_{\text{SB}}^{\ast }/\sum _S^{\mathrm{2}}+\widehat{U}{T}_{\text{UB}}^{\ast }/\sum _U^{\mathrm{2}}+\widehat{V}{T}_{\text{VB}}^{\ast }/\sum _V^{\mathrm{2}},\\ {\displaystyle }& {\displaystyle }{Y}_{\mathrm{c}}=\widehat{S}{T}_{\text{SC}}^{\ast }/\sum _S^{\mathrm{2}}+\widehat{U}{T}_{\text{UC}}^{\ast }/\sum _U^{\mathrm{2}}+\widehat{V}{T}_{\text{VC}}^{\ast }/\sum _V^{\mathrm{2}}.\end{array}$$

The least-squares solution is then

$$\begin{array}{ll}{\displaystyle }\left[\begin{array}{c}\widehat{B}\\ \widehat{C}\end{array}\right]& {\displaystyle }={\left[\begin{array}{cc}L& K\\ K^{\ast }& M\end{array}\right]}^{-\mathrm{1}}\left[\begin{array}{c}{Y}_{\mathrm{b}}\\ Y_{\mathrm{c}}\end{array}\right]\\ \text{(C3)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{LM-K{K}^{\ast }}}\left[\begin{array}{cc}M& -K\\ -K^{\ast }& L\end{array}\right]\left[\begin{array}{c}{Y}_{\mathrm{b}}\\ Y_{\mathrm{c}}\end{array}\right].\end{array}$$

This inversion is problematic where $LM-K{K}^{\ast }$ is small or zero, which is the case for small wavelengths (when k and l are large) or for topography or slipperiness perturbations which are aligned in the direction of ice flow (k=0). Short-wavelength bed features and features aligned with ice flow are problematic because they cause flow disturbances in the ice which do not reach the surface in a measurable way. They can therefore not be inverted from the surface data.

To avoid this problem with an ill-conditioned inverse, Thorsteinsson et al. (2003) use a truncated version of $L^{\prime }{M}^{\prime }-K{K}^{\ast }$ as a filter, F, to remove problematic wavelengths. Their filter is

$$F=\left\{\begin{array}{ll}(LM-K{K}^{\ast })/P& \text{if}(LM-K{K}^{\ast })\le P,\\ \mathrm{1}& \text{if}(LM-K{K}^{\ast })>P,\end{array}\right.$$

where

$$\begin{array}{}\text{(C4)}& P=\text{max}\mathit{\left\{}\right(|LM-K{K}^{\ast }|\left){\stackrel{\mathrm{‾}}{C}}^{\left(p_{\text{filt}}\right)}\mathit{\right\}}\end{array}$$

and p_filt≤0. This filter allows through all wavelengths where $LM-K{K}^{\ast }$ is larger than P but gradually filters out other smaller wavelengths. Smaller values of P give more detail but may over-fit the surface data due to errors. Larger values of P under-fit the data and may leave out features actually represented by the data.

The filtered least-squares solution is then

$$\begin{array}{}\text{(C5)}& \left[\begin{array}{c}\widehat{B}\\ \widehat{C}\end{array}\right]={\displaystyle \frac{F}{LM-K{K}^{\ast }}}\left[\begin{array}{cc}M& -K\\ -K^{\ast }& L\end{array}\right]\left[\begin{array}{c}{Y}_{\mathrm{b}}\\ Y_{\mathrm{c}}\end{array}\right].\end{array}$$