the Creative Commons Attribution 3.0 License. The Cryosphere Discussions

An optimal estimation method for simultaneously determining both basal slipperiness and basal topography from variations in surface flow velocity and topography along a flow line on ice streams and ice sheets is presented. We use Bayesian inference to update prior statistical estimates for basal topography and slipperiness using surface measurements along a flow line. Our main focus here is on how errors and spacing of surface data affect estimates of basal quantities and on possibly aliasing/mixing between basal slipperiness and basal topography. We find that the effects of spatial variations in basal topography and basal slipperiness on surface data can be accurately separated from each other, and mixing in retrieval does not pose a serious problem. For realistic surface data errors and density, small-amplitude perturbations in basal slipperiness can only be resolved for wavelengths larger than about 50 times the mean ice thickness. Bedrock topography is well resolved down to horizontal scale equal to about one ice thickness. Estimates of basal slipperiness are not significantly improved by accurate prior estimates of basal topography. However, retrieval of basal slipperiness is found to be highly sensitive to unmodelled errors in basal topography.


Introduction
Indirect inference of basal shear stress distributions of ice streams from measurements of surface data has been used successfully in the past to study the role of basal conditions in controlling the flow of ice streams (MacAyeal, 1992(MacAyeal, , 1993MacAyeal et al., 1995;Vieli and Payne, 2003; Correspondence to: G. H. Gudmundsson (ghg@bas.ac.uk) 2004Gudmundsson (ghg@bas.ac.uk) , 2006. Simultaneous retrieval of both basal topography and basal lubrication was done by Thorsteinsson et al. (2002), and a Bayesian framework for surface-to-bed inversion was outlined by Gudmundsson (2004). To date, limited attention has been paid to some more general aspects of surface inversion on glaciers. For example it remains unclear how the accuracy of estimates of basal quantities is affected by the spatial distribution of surface measurements and measurement errors. Possible mixing effects, such as the effects of basal topography on retrieved slipperiness distributions have also not received much attention.
The aim of this paper is to clarify the effects of data density and data quality on indirect estimates of basal topography and basal slipperiness on ice streams. We do this by considering three different aspects of the retrieval, 1) errors due to mixing between basal topography and basal slipperiness in retrieval and mixing between different frequency components (mixing errors), 2) the spatial and frequency resolution of the retrieval, and 3) the number of resolved basal quantities per measurement site at the surface.
We use a Bayesian approach to this inverse problem and determine the bedrock profile (b(x)) and the spatial slipperiness distribution (c(x)) which maximise the conditional probability P (b, c|s, u, w), where s(x), u(x), and w(x) are the surface topography, and the horizontal and vertical surface velocity component, respectively. The forward model is based on perturbation solutions to the full Stokes system (Gudmundsson, 2003). Consequently our results are only strictly valid in the limiting case of linearised flow perturbations. Numerical estimates of the strengths of non-linear effects in the forward model can be found in Raymond and Gudmundsson (2005).
G. H. Gudmundsson and M. Raymond: Surface inversion on ice streams The method used here differs in a number of aspects from previous work done by other authors on the surface-to-bed inversion problem on glaciers and ice sheets. We formally introduce prior information about basal properties and use the surface measurements in combination with the prior information to give an updated estimate of basal properties. This new estimate of basal properties is given in terms of a Gaussian probability density function defined through a maximum a posteriori (MAP) solution and its covariance. In the sense that the solution maximises the conditional probability P (b, c|s, u, w), it can be regarded as an optimal estimate of basal properties. A further example for optimal estimation in ice sheet modelling can be found in (Arthern and Hindmarsh, 2003).
Although not always stated explicitly, some sort of a priori expectation about "sensible" basal properties has been used in most, if not all, previous work. For example the particular minimising norm used by Truffer (2004) reflects the reasonable expectation that basal properties are in some ways spatially correlated and that surface data can only be expected to give information about basal properties over restricted set of spatial scales.
As pointed out by (MacAyeal et al., 1995), for most inverse problems direct minimisation of the misfit between measured and modelled data does not always lead to a satisfactory solution. In the control method introduced by MacAyeal (1993) and used by number of other authors (e.g., Vieli and Payne, 2003;Joughin et al., 2004;Larour et al., 2005), this problem is avoided by starting with an initial guess judged to be reasonable by the modeller and by not fully minimising the misfit between measured and modelled data. Not fully minimising the misfit is a way of ensuring that the solution norm does not grow beyond bounds. In effect this limits the set of admissible solutions and can be thought of as a form of a constraint on the solution.
Although hitherto not commonly used in glaciology, the Bayesian approach to inverse problems applied here is often used to solve geophysical inverse problems, and we refer the reader to books by Menke (1989), Rodgers (2000), and Tarantola (2005) for a general description. Bayesian inversion is a statistical inference method where indirect measurements of some quantity are used in combination with earlier estimates of that quantity to arrive at a new and improved estimate. The estimate is optimal in the sense that it maximises the conditional probability P (x|y) where x is the quantity to be estimated, and y are indirect measurements of x. The indirect measurements y of x are related to x through y=f (x) where f is referred to as the forward model. Here, in our particular case, the indirect measurements (y) are measurements of velocity and topography along the upper surface of a glacier, the quantity x to be estimated is basal topography and basal slipperiness.

Theoretical framework
We consider the problem of determining basal slipperiness and basal topography from measurements of surface topography and surface velocity on glaciers. The forward problem consists in solving the full set of the linear momentum equations σ ij,j +f i =0, the angular momentum equations σ ij −σ j i =0, and the incompressibility condition v i,i =0 for a linear media where˙ ij =τ ij /(2η) and for a linear sliding law u b =cτ b . In these equations σ ij are the elements of the Cauchy stress tensor, f i the elements of the body force, v i the velocity components,˙ ij are strain rates, τ ij are deviatoric stresses, η the ice viscosity, u b the basal sliding velocity, τ b the basal shear stress, and c the basal slipperiness. In addition to these equations we have the kinematic boundary conditions at the upper and lower boundaries. We consider steadystate conditions and ignore accumulation and ablation.
The model domain is that of a slightly perturbed uniformly inclined slap with periodic boundary conditions. The perturbations that we consider here are perturbations in basal topography and basal slipperiness. Unless otherwise specified, the amplitude of a basal topography perturbation is given as a fraction of the mean ice thickness. Correspondingly, the amplitude of a basal slipperiness perturbation is given as a fraction of the mean basal slipperiness. For example, a (fractional) basal slipperiness perturbation with an amplitude equal to 0.1 corresponds to a 0.1c perturbation in basal slipperiness, wherec is the mean basal slipperiness. We havē u b =cτ b whereū b andτ b are the mean basal sliding velocity and the mean basal shear stress, respectively.

Notation
Measurements are available at discrete points, and we refer to the set of all surface measurement values as the measurement vector. When referring to this (finite) set of surface measurements, the corresponding vectors are written in bold face.
The variables to be estimated are the vectors b and c giving the bed topography and basal slipperiness along the profile. We lump these vectors together into one system state vector x and, think of x as consisting of two subvectors x 1 and x 2 where x 1 =b and x 2 =c. Hence x=(b, c) H , where H denotes the Hermitian transpose. We use the Hermitian transpose because the inversion is done in Fourier space where all variables are complex.
Surface data is lumped together into the measurement vector y. The set of surface data considered here are measurements of surface topography (s), and the horizontal and vertical velocity components u and w, respectively. For notational convenience we split y into three subvectors and write y 1 =s, y 2 =u, y 3 =w. Hence y=(s, u, w) H . We also write The set of all parameters entering the problem that are not solved for as a part of the inversion procedure is denoted by p. These parameters include, for example, those related to ice rheology.
A tilde above a quantity denotes an a priori estimate while a hat above a quantity denotes an a posteriori estimate. Thus x is the true system state whilex, andx are the a priori and the a posteriori estimates, respectively. Measurement values are denoted by a hat above the corresponding quantity. We write for exampleŝ=s+ s whereŝ is the surface measurement vector, s is the true surface and s the error in the surface measurement.

Linear forward problem
The relationship between surface measurements and basal variables is written as where f is the forward model. We have assumed here that the forward model is perfect in the sense that it correctly describes the relationship between basal (x) and surface (y) fields when exact estimates of model parameters (p) and basal quantities are available.
In Fourier space the relationship between horizontal and vertical surface velocity components (u and w), surface topography (s) and the basal variables b and c can be described through transfer functions. Here we use the transfer functions given by Gudmundsson (2003). Surface and bedrock variables are related through where the summation convention is used, and y i T x j are transfer functions. It is useful to use extended transfer functions, defined as These matrices have the dimensions m×2n where n is the number of locations along the bedline at which c and b is to be estimated, and m is the number of measurement sites at the surface where s, u and w have been measured. In general, no notational difference is made between the same function in spatial and frequency domain. In the discrete case, the forward Fourier transform of the corresponding measurement vector s is given by Fs where F is the unitary discrete Fourier transform matrix. The discrete forward Fourier transform of the covariance matrix C is given by FCF H where H denotes the Hermitian transpose.

Description of errors
Covariance matrices are denoted by the upper case letter C. The covariance matrices considered here will almost all be error covariance matrices, where the expectation value operator is taken about the true value of some quantity. The error in surface topography measurements is, for example, described through the covariance matrix or simply as We assume that measurement errors are normally distributed and uncorrelated. The error in the prior estimate is described by the covariance of the ensemble of states about the prior state We assume that the prior estimates of bed topography and basal slipperiness are mutually uncorrelated and write C xx on the form We know of no studies were the spatial correlation in basal topography or basal slipperiness have been estimated. It seems, however, unlikely that basal topography and basal slipperiness are spatially uncorrelated. A reasonable model for the prior basal topography covariance is a stationary first-order auto-regressive process (Box et al., 1994) where the prior estimate of the basal topography at location i is related to that at location i−1 through where b is uncorrelated with variance σ 2 . It follows (Box et al., 1994) that the elements of the covariance matrix C bb are where the decay length scale λb is related to φb through and the variance is The two parameters λb and σ 2 b define the covariance model. We refer to σb as the error in prior basal topography. In what follows this error will usually be defined in terms of the mean ice thicknessh. For example, for a 10 % error in prior basal topography we have σ b =0.1h. Do to the to lack of corresponding studies it is difficult to give any well justified estimates for λb. Over spatial scales where the surface data contains enough information about basal properties for them www.the-cryosphere.net/2/167/2008/ The Cryosphere, 2, 167-178, 2008  The blue lines are the retrieved profiles using Eq. (17). Lines in red show the total retrieval error estimate defined as the square root of the diagonal elements of the C xx covariance matrix (Eq. 18). The lines in cyan show the a priori retrieval error (see Eq. 24). the be well resolved, the value of λb has little or no effect on the retrieval. Over spatial scales where the surface data does not lead to any significant improvements in the prior estimate of basal properties, the value of λb, on the other hand, significantly affects the spatial correlation ("smoothness") of the retrieval. With the exception of the results shown in Fig. 1, which is an illustration of one particular synthetic inversion experiment, all the results given below are not affected by the exact value of λb An identical model is used for the prior basal slipperiness covariance as for the basal topography, and the corresponding model parameters are denoted by σc and λc. A 10 % error in basal slipperiness corresponds to σc=0.1c.
One of the reasons for using this model for the prior covariance C xx is that analytical expression for its inverse (C −1 xx ) and its Fourier transform (FC xx F H ) can easily be determined (Box et al., 1994). C xx is the only full covariance matrix entering as an input variable to the problem, and having analytical expressions for its inverse, and its Fourier transform, significantly speeds up the inversion procedure. The retrieval error is described through the covariance matrix In general we expect C xx to be a full covariance matrix and the estimates of basal topography and basal slipperiness to be correlated.

Inverse problem
The retrieval (x) is a function of surface data vectors (ŝ,û, w), the prior estimate of the system state (x), and the set of estimated model parameters (p), that iŝ We seek to determine a system state that maximises the conditional probability P (b, c|u, v, w). Using Bayes' theorem this probability can be written as is the prior estimate of basal properties. Standard arguments show that the maximum a posteriori (MAP) estimate of the system state is given by

Error and sensitivity analysis
Linearising Eq. (15) shows that the difference between the estimate of the system state (x) and the true system state (x) is given bŷ wherê is the sensitivity of the retrieval (x) to the true system state (x), and is the sensitivity of the retrieval to forward model parameters, and is the sensitivity of the retrieval to surface data with y 1 =s, y 2 =u, y 3 =w.
The first term on the right-hand side of expression Eq. (20) is the retrieval error due to inaccuracies in the prior. The corresponding covariance matrix is We refer to this error component as the "prior retrieval error", r .
The second term on the right-hand side of expression Eq. (20) is the retrieval error due to errors in the forward model parameters and the sensitivity of the retrieval to the forward model parameters. We refer to this error component as the "forward-model parameter retrieval error", f . Here we assume that errors in forward model parameters are small. If not, the corresponding model parameters should be included as a part of the system state.
The third term on the right-hand side of Eq. (20) is the retrieval error due to measurement errors ( n ) and the sensitivity of the retrieval to measurements. The (total) retrieval error given by Eq. (18) is the sum of the retrieval error due to measurements errors and the prior retrieval error.
For the MAP estimate given above and the definition of x A x , it follows, using the chain rule and the assumption that measurements are unbiased and errors uncorrelated, that We refer toxA x as the averaging kernel matrix. The matrix x A x is, in general, not a symmetrical matrix. It is a 2n×2n matrix, where 2n is the number of elements in x.xA x can be written aŝ where each of the blocksbA b ,bA c ,ĉA b , andĉA c are n×n. ThebA c matrix gives the sensitivity of the bedrock retrieval to the (true) basal slipperiness distribution, while theĉA b matrix gives the sensitivity of the estimate of basal slipperiness to the (true) bedrock topography. These two averaging kernel matrices will be referred to as mixing kernel matrices. Clearly, for a good retrieval these elements of the mixing kernel matrices should be as small as possible. The elements ofb A b give the sensitivity of the bedrock-topography retrieval to those of the true bedrock topography. SimilarlyĉA c is the sensitivity of the slipperiness retrieval to the actual basal slipperiness distribution. We refer to thebA b and theĉA c kernels as the direct averaging kernel matrices. For good retrieval there should be as little mixing as possible mixing between frequencies while the sensitivity of corresponding frequencies in the true and estimated system state should be as close as possible to one, that is ideally the direct kernelsbA b and c A c should be unit matrices.
The sensitivity of the retrievalx to the prior estimatex is defined aŝ From Eq. (17) we find and using either Eq. (20) to calculatexAx together with Eq. (25), or using the above equation in combination with Eqs. (18) and (25), it follows that giving the relationship between the sensitivity of the retrieval to the system state and the prior estimate, and showing that the retrieval is sensitive to either the true system state or the prior estimate. Hence, deviations of the averaging kernel ma-trixxA x from the unit matrix is a measure of how sensitive the retrieval is to the prior estimate as opposed to the surface data. Taking the trace of Eq. (29) we can define the number of model parameters resolved by the measurements as trxA x and the number of model parameters resolved by the prior information as trxAx. The number of basal quantities resolved by the surface measurements (d s ) is therefore i.e. the trace of the averaging kernel matrix.

Results
As an illustration of the method we start by showing results of an inversion using synthetic data. We stress that apart from this one example, none of the following results depend on the use of synthetic data. The forward calculation was done for an ice stream with a mean surface slope of 0.005, and a slip ratio of 500, where the slip ratio is defined as the ratio between mean basal sliding velocity and forward velocity due to shearing through the ice column. Spatially uncorrelated and normally distributed data errors were added. Surface topography errors were 1% of mean thickness, and errors in both horizontal and vertical velocity components were set at 1% of mean surface velocity. The spacing between measurements was equal to one mean www.the-cryosphere.net/2/167/2008/ The Cryosphere, 2, 167-178, 2008 ice thickness. The prior for the basal topography and the basal slipperiness was modelled as a covariance-stationary first order auto-regressive model with zero mean as described in the previous section. Both decay lengths (λb and λc) were set at 25h (see Eq. 11). The error was 10% of the respective mean values, i.e. σb=0.1h and σc=0.1c (see Eq. 12). Synthetic prior distributions (x) were generated usingx=C where is a sequence of pseudo random numbers with mean zero and unit variance. Gaussian perturbations in true basal topography and true basal slipperiness were prescribed. Non-dimensional variables are used. All lenghts are normalized with the mean ice thickness. Stresses are normalized with the mean basal shear stress, and velocities by mean deformational velocity. These scaling follow Gudmundsson (2003). As before the ratio between mean sliding velocity and basal shear stress is referred to as the basal slipperiness. In non-dimensional variables this ratio is also equal to the ratio between mean sliding velocity and mean deformational velocity, that is, to the slip ratio. Hence, in the following the terms "slip ratio" and "basal slipperiness" are interchangeable.
An inversion of synthetically generated surface data with added errors was performed using Eq. (17). As Fig. 1 shows, the retrieved bedrock perturbation (blue line, Fig. 1a) traces the true perturbation (black line) quite accurately. The retrieval errors are given by Eq. (18). These errors are in general correlated and for that reason somewhat difficult to visualise. In Fig. 1, the square root of the diagonal elements of C xx is used to bracket the range of possible retrievals (red lines). As explained above this error can be thought of as the sum of the prior retrieval error ( r ) given by Eq. (24) and retrieval error due to measurement errors ( n ). For the basal topography retrieval, the prior retrieval error component, shown as cyan line in Fig. 1a, is a small fraction of the total retrieval error (red line, Fig. 1a). Hence, the error in basal topography retrieval is mostly due to surface measurement errors rather than errors in prior estimates of basal properties.
In comparison to the estimate of basal topography, the retrieved basal slipperiness perturbation (blue line, Fig. 1b) is a much poorer estimate of the true basal perturbation (black line, Fig. 1b). As Fig. 1b shows, the main contribution to the total retrieval error (red line) is due to errors in the prior estimate (cyan line). It follows that increasing the accuracy of the surface data somewhat will not lead to significant improvements in slipperiness retrieval and that the main limitation to further improvement in slipperiness retrieval is insensitivity of the retrieval to surface data.

Spatial resolution and mixing effects
Inversion of synthetic data of the type shown in Fig. 1 can be instructive and helpful in identifying problems with the inversion method. However the usefulness of such an inver-sion exercise is limited by the need to prescribe a particular form of basal perturbation.
A general description of the sensitivity of the retrieval to the system state is given by the averaging kernel matrixxA x (see Eqs. 21 and 25). Considering the averaging kernel matrix in frequency space we note that for each of the blockŝ b A b ,bA c ,ĉA b , andĉA c , the i-th row represents the sensitivity of the retrieval to a system state consisting of white phaseless noise. If there is no mixing between frequencies, that is if one frequency in the retrieval is only related to that same frequency of the system state, then the only non-zero elements of the i-th row of these blocks is at location i of the corresponding block. It follows that for each of the blockŝ b A b ,bA c ,ĉA b , andĉA c , the numerical difference between the element at location i, and the sum over all elements of the i-th row, can be taken as a measure of frequency mixing. Similarly, for each of the blocksbA b ,bA c ,ĉA b , andĉA c the i-th columns give the frequency spread, i.e. the sensitivity of all frequencies of the retrieval to one single frequency of the system state.
A further undesirable mixing effect is represented by the blocksbA c and theĉA b . These give, respectively, the sensitivity of the basal topography retrieval to true basal slipperiness, and the sensitivity of the retrieved basal slipperiness to true basal topography. For a good retrieval these matrix elements should not only be small in comparison to unity, they should also be small in comparison to the diagonal elements along the same rows of the averaging kernel matricesbA b andĉA c . Figure 2 shows the sensitivity of the retrieval to true system state in frequency space. The curves in the figure follow from Eqs. (25) and definition (26). The difference between solid and dashed lines of same colour represent frequency mixing effects. The green and the cyan lines represent mixing between basal topography and basal slipperiness in the retrieval. The figure was calculated assuming errors in surface measurements of topography and velocity equal to 1% of mean values of ice thickness and surface speed. For an ice stream with thickness of 1000 m and surface speed of 1 m d −1 this translates to 10 m elevation errors and 0.01m d −1 errors in horizontal and vertical velocity components. Airborne measurements of surface topography can easily reach this kind of accuracy (Vaughan et al., 2006), and both repeated annual surveying of stakes positions using GPS techniques and InSAR velocity measurements give surface velocities within or comparable to this level of error (Joughin, 1995;King, 2004;Dach et al., 2009).
The dashed-lined red curve in Fig. 2 shows the size of the diagonal elements of thebA b averaging kernel matrix. These are close to unity showing that bedrock topography is well resolved down the wavelengths comparable to one mean ice thickness. The solid red curve is the sensitivity of the retrieval at given wavelength to white noise. The small show the sensitivity of the retrieval as a function of wavelength to all wavelengths of the true basal conditions. The difference between the solid and the dashed lines of same colour is a measure for undesirable frequency mixing effects. The red lines show the sensitivity of the estimate of basal topography to the true basal topography, and the blue lines the sensitivity of estimated basal slipperiness to true basal slipperiness. The green and the cyan lines represent further undesirable mixing effects where bedrock perturbations affect estimates of basal slipperiness (cyan lines) or where estimates of bedrock topography are contaminated by variations in basal slipperiness. The figure was calculated using a surface slope of 0.005 and mean slip ratio of 500. Surface measurement locations were one ice thickness apart and surface data errors were 1% of ice thickness and mean surface velocity. prior values for basal topography and basal slipperiness were modelled as first-order auto-regressive processes with decay length scales of 25h and 10% errors. difference between the solid and the dashed red curves implies that frequency mixing is not a significant problem for bedrock topography retrieval.
The small difference between the dashed and the solid blue lines in Fig. 2 shows that frequency mixing does not pose a problem for basal slipperiness retrieval either. However the spatial resolution is somewhat limited with sensitivity greater than 0.5 only for wavelengths larger than about 40-50 ice thicknesses.
Mixing between different types of basal perturbations is represented by the green and cyan curves in Fig. 2. These dashed lines show the main diagonal elements of thebA c and theĉA b averaging kernels, respectively, while the solid lines are the sums of the rows. Significant mixing (>0.2) is only found in theĉA b kernel (solid cyan line in Fig. 2) and is limited to wavelengths within a fairly narrow range of 20 to 40 ice thicknesses. As can be seen from the substantial differ- The upper panel shows in blue the sensitivity of basal topography (b) retrieval to point perturbation in basal topography located at x=0, and in green the sensitivity of basal slipperiness (c) retrieval to point perturbation in slipperiness. The lower panel shows the degree of mixing between basal topography and basal slipperiness in the retrieval. The blue curve in the lower panel is the sensitivity of the basal topography retrieval to basal slipperiness, and the green curve is the sensitivity of the basal slipperiness retrieval to point perturbation in basal topography. Values of all model parameters are same as those used in Fig. 2. ence between the dashed and the solid cyan lines, this mixing between basal slipperiness retrieval and basal topography is combined with frequency mixing. Hence, sinusoidal perturbations in basal topography are aliased into basal slipperiness perturbations (b to c mixing) at different wavelengths (frequency mixing).
From Fig. 2 it is not clear which frequencies of the system state contribute to frequency mixing. It is also difficult to visualise the spatial resolving power of the retrieval method from a frequency-space representation such as the one in Fig. 2. A different way of looking at the averaging kernels is to consider the "point spread functions" shown in Fig. 3. The point spread functions represent the spatial, as opposed to the frequency, sensitivity of the retrieval to a localised point perturbations in system state. Figure 3 shows four point spread functions corresponding to the sensitivity of the basal topography retrieval to point perturbations in basal topography (blue curve labelled bb in Fig. 3a) and basal slipperiness (blue curve labelled bc in Fig. 3b), and the sensitivity of the basal slipperiness retrieval to point perturbations in basal slipperiness (green curve labelled cc in Fig. 3a) and basal topography (green curve labelled cb in Fig. 3b). Values of all relevant parameters such as average slope, data errors and data spacing are the same as in Fig. 2. The strong differences in spatial resolution of the basal topography retrieval in comparison to that of the basal slipperiness retrieval is clearly seen in Fig. 3a. As the figure shows, a point perturbation in the system state of b is retrieved almost perfectly while the retrieval of a point perturbation in basal slipperiness retrieval is broad.
The degree of mixing between different types of basal perturbations is shown in Fig. 3b. As Fig. 2 showed in a different way, this type of mixing is not strong. In particular there is almost no mixing with the basal slipperiness in the basal topography retrieval. The green curve in Fig. 3b shows that point perturbation in b leads to some reduction in slipperiness in upstream direction from a bedrock disturbance and a similar increases in downstream direction.

Number of resolved basal quantities
Our discussion above on spatial resolution and mixing effects was limited to one particular set of parameters defining data errors and data spacing. All surface data errors were 1% of mean values or ice thickness and horizontal speed, and prior errors were set at 10% of mean thickness and mean slipperiness with decay length scale equal to 25 ice thicknesses. Furthermore, the spacing between measurement sites at the surface was equal to one mean ice thickness. We will now investigate the effects of surface data errors and spacing on the number of resolved basal quantities as given by Eq. (30).
A convenient measure of information on basal properties gained by surface measurements is the number of resolved basal parameters per measuring site at the surface. As Eq. (30) shows the total number of basal quantities resolved by surface measurements equals the trace of thexA x averaging kernel matrix. For both basal topography and basal slipperiness the maximum number of resolved parameters is n, and the maximum total number of resolved parameters therefore 2n. The number of measuring sites at the surface is m, and at each of these sites we have three data values giving the surface elevation and the horizontal and vertical velocity component at that site. In the following we assume that the number of measurement sites equals the number of locations along the bedline where estimates of basal quantities are to be updated, that is m=n. The number of resolved basal topography (d b s ) and basal slipperiness (d c s ) quantities per measurement site is then and respectively. If, for example, the number of resolved basal topography quantities per measurement site is close to unity, all of the surface data sites are contributing significantly more to the estimate of basal topography than the prior information. If this number is equal to 0.5, the sensitivity of the retrieval to surface data is equal to its sensitivity to prior information.
It should be noted that there are number of other possible ways of quantifying the information content of the surface measurements. One could for example define the information content of surface measurements in terms of the reduction in entropy (Shannon and Weaver, 1949;Togneri and DeSilva, 2003), in which case the difference between the covariance of the MAP solution and the covariance of the prior is taken as a measure of how much the measurements have affected our estimate of the system state. We calculated the Shannon information content but found, as have other authors (e.g. Purser and Huang, 1993), that this measure of information has some undesirable properties. For example, despite being only of local influence, a single precise measurement greatly alters this measure of information. Figure 4 shows how data errors affect the number of resolved basal quantities per surface measurement site as a function of distance between sites. Data errors are shown as percentages of corresponding mean values. A 0.1% data error implies, as an example, a 0.001h error in surface topography and a 0.001ū error in both horizontal and vertical velocity components, whereū is the mean surface speed. Ideally, the number of resolved quantities per site would be close to unity. Note also that a slope of zero in Fig. 4 represents a situation where the effect of reducing the spacing between measurement sites is a directly proportional increase in the total number of resolved basal quantities. As Fig. 4 shows the number of resolved basal topography quantities per site only starts to drop significantly down from unity with reducing site spacing for spacing corresponding to about one ice thicknesses, with exact values depending on the size of data errors (Fig. 4, red curves). For basal slipperiness, on the other hand,d c s is only larger than about 0.7 for spacing larger than 10h if data errors are less than 0.1% (Fig. 4, blue curves). Figure 5 showd b s andd c s when one of the surface data types (s, u, or w) is not known. Of the three data sets s, u, and w, not having information about the surface topography (s) has the largest detrimental effect on the b retrieval. However, as the figure demonstrates knowledge of the surface topography is far from critical for the retrieval of basal topography. Almost identically many basal topography parameters can be resolved in the absence of information about surface topography from knowing the surface velocity vector alone. It is only for distances above about 20 ice thicknesses that not knowing the surface topography starts to significantly degrade the b retrieval. A further interesting aspect of Fig. 5 is how unimportant the vertical velocity component (w) is for the c retrieval. In fact, not knowing w has such a small effect ond c s that the corresponding blue coloured dash-dotted curve in Fig. 5, showingd c s for s and u but not w known, is not visible under the thicker solid blue curve givingd c s for all of s, u and w known.  Fig. 6. Number of resolved basal topography (red curves) and basal slipperiness (blue curves) quantities as a function of distance between measurement sites at the surface for different error estimates for the prior basal topography. The thick and the thin curves were calculated given 1% and 30% uncertainty, respectively, in the prior basal topography.
The basal slipperiness retrieval can be expected to be improved when accurate prior information on basal topography is available. Figure 6 shows the number or basal quantities resolved by surface measurements for two different error estimates of prior basal topography. Comparison of the two blue lines in that figure shows that improving the estimates of prior basal topography from 30% error level to 1% error does indeed increase the number of resolved basal slipperiness quantities. However, that improvement is modest and the increase ind c s no larger than about 0.1 (Fig. 6). We conclude that a good prior knowledge of basal topography is not an essential requirement for successful retrieval of basal slipperiness.
The number of basal topography quantities resolved by the prior clearly must increase with decreasing prior errors. Hence the significant drop ind b s values with reduced errors in prior basal topography seen in Fig. 6. Interestingly, the sensitivity of the basal topography retrieval to surface data is still around 0.2 to 0.3 for site distances above one mean ice thickness for only 1% errors in prior basal topography (see Fig. 6). Hence, it is important to solve for both b and c even when accurate independent estimates of b are available. Not solving for b when basal topography is known corresponds to artificially suppressing the sensitivity of the b retrieval to surface data to zero, whereas in fact, as Fig. 6 shows, the actual sensitivity for just 1% prior errors in b is around 0.2 to 0.3. Ignoring this sensitivity will inevitably cause a bias in the c estimate. The magnitude of this bias can be estimated by considering the mixing terms of the kernel matrixxAx. From Eq. (29) and using Eq. (26) we find that sensitivity of the c retrieval to prior estimates of b is the negative ofĉA b . As we have seen this mixing term can be as large as about 0.3 for wavelengths between around 20 to 50 ice thicknesses with most of that contribution due to frequency mixing (see solid cyan curve in Fig. 2). We note that Fig. 2 was calculated for a σb=0.1h and that the exact errors in c retrieval due to incorrect error estimates for the prior basal topography will change somewhat as a function of σb. Inspection of x A x and its blocks shows that with decreasing prior errors in basal topography the sensitivity of the c retrieval to perturbations in true basal topography increases sharply and can easily become much larger than unity.
The sensitivity of the retrieval of basal slipperiness and basal stress to errors in basal topography was investigated by Joughin et al. (2004). They solved for basal shear stress only, without allowing for any variations in basal topography. Contrary to our findings they concluded that errors in basal shear stress estimates did not increase markedly with errors in basal topography. Similar conclusions were drawn by Vieli and Payne (2003) who inverted for basal slipperiness with, and without, a simultaneous inversion for bedrock topography, and found that in both cases the resulting basal slipperiness was quite similar. It is possible that the root of some of this difference might be in the use of different forward models. Joughin et al. (2004) solved the shallow ice stream equation for non-linear media while we solve the full Stokes system for linear media. However, it is also possible that this difference is related to difference in the way the sensitivity experiments were performed. Joughin et al. (2004) solved the diagnostic equation for each change in input parameters. The prognostic equation was not used to ensure that both the rate of surface elevation changes and the surface mass balance were within set limits. This is an approach followed by number of other authors (e.g., Vieli and Payne, 2003). In this study we have assumed zero surface mass balance and zero rate of surface change with time and have accordingly calculated the steady state surface topography and surface velocities corresponding to zero surface mass balance for any given change in basal quantities. It seems likely that the sensitivities of retrieved basal quantities to errors in basal topography when a) surface topography is held constant, and when b) corresponding steady state surface topography is calculated, are not the same. The differences between our conclusions and those of Joughin et al. (2004) could, hence be due to differences in the way the sensitivity of retrieval to errors in basal topography is defined.

Limitations and possible extensions
A key feature of the retrieval method is the use of analytical transfer functions, that describe the effects of smallamplitude perturbations in basal conditions on surface geometry and surface velocity, as a forward model. The transfer functions used here are solutions to the full Stokes system and only valid for Newtonian media and a linear sliding law. However, the effects of a non-linear sliding law on surface-to-bed retrieval on ice streams, over horizontal spatial scales that are large compared to the mean ice thickness, could be studied by using the transfer functions given in Gudmundsson (2008). Doing so would require no essential modifications to the method as presented here. Including finite-amplitude effects in combination with non-linear ice rheology and non-linear sliding law requires a fully numerical treatment. This non-linear retrieval problem can, for example, be solved using an iterative procedure where forward model derivatives are approximated by the analytical transfer functions. The iterative step of such a method would be almost identical to the method presented and used here. This extension of the method to non-linear rheology, non-linear sliding law, and finite-amplitude perturbations is the subject future paper (Raymond and Gudmundsson, 2008).
Transfer functions for non-linear rheology and non-linear sliding have been calculated numerically by Raymond and Gudmundsson (2005). Qualitatively, the effects of non-linear rheology and non-linear sliding on the retrieval can be estimated by comparing transfer functions calculated for both the linear and non-linear situations. For high slip ratios (>5), the effect of increasing the value of the stress exponent in Glen's flow law is an increase in the ratio between surface and bed topography (see Fig. 8b in Raymond and Gudmundsson, 2005, for more detail). For the ratio between surface topography and basal slipperiness the effect is, over most wavelengths, the opposite and the (non-dimensional) ratio between perturbations in surface topography and fractional slipperiness decreases with increasing value of the stress exponent (see Fig. 9b in Raymond and Gudmundsson, 2005, for more detail). This suggests that non-linear rheology may improve the retrieval of basal topography but degrade slipperiness retrieval. However, the picture is complicated by the fact that changing the value of the stress exponent also affects surface velocity response to basal perturbations, and a definitive answer must await a fully non-linear treatment.
The method is here only used for data along a flow line and we have not considered the consequences of the transverse dimension on retrieval. In particular, we empathise that the bounds given above on the spatial resolution of basal slipperiness variations have been obtained for a flow-line inversion where transverse effects have been ignored. Perturbations in basal slipperiness of final transverse width give rise to patterns of horizontal divergence and convergence (see Figs. 2, 5c and 5d in Gudmundsson, 2003). It is possible that in many situations the resulting disturbance in horizontal flow velocities are large enough for perturbations in basal slipperiness of finite width to be extracted from surface data with higher spatial resolution than suggested by the one-dimensional treatment given here. Based on the transfer characteristics (see Fig. 5b in Gudmundsson, 2003) it seems, for example, quite possible that a transversal sinusoidal perturbation in basal slipperiness, i.e. with a crest that is aligned with the flow, can be resolved from surface data down to a (transverse) wavelength considerably shorter than the limit of about 50h given here for longitudinal perturbations. Further support for this expectation comes from Schoof (2004) who found that a narrow transversal perturbation in basal yield-stress on a perfectly plastic ice stream can significantly affect the horizontal velocity component at the surface (see Fig. 3a in Schoof, 2004).
Finally, we have not considered the effects of finiteamplitude perturbations in basal properties, or the effects of form drag on the transmission of basal disturbances to the surface. Our purpose here has been to give as general a description of the surface-to-bed inverse problem as possible using an analytical approach. Finite-amplitude perturbations are inherently non-linear, the superposition principle does not apply and making general statements based on an analytical approach is difficult. However, we do not expect finite-amplitude effects, or the effects of form drag, to change our results. Raymond and Gudmundsson (2005) found finite-amplitude effects not to significantly affect the bed-to-surface transmission characteristics for small to moderate basal amplitudes. For example, for sinusoidal basal perturbations with a wavelength λ=5h, the analytical solutions are correct to within a few percent for normalised bedrock amplitudes ( b/h) and fractional slipperiness ( C/C) amplitudes up to 0.3h for a slip ratio of 1000 (see Fig. 3 in Raymond and Gudmundsson, 2005). For large wavelengths the agreement is even better. These findings are in sharp contrast to conclusions drawn by Schoof (2002) who argued that some of the assumptions behind the analytical approach in (Gudmundsson, 2003) start to be violated for bedrock amplitudes as small as 0.2 m for 1000 m thick icestreams.

Summary and conclusions
We have shown how basal topography and basal slipperiness can be retrieved simultaneously from surface measurements without any significant frequency mixing or mixing between these two different types of basal perturbations. No smoothing of surface data is needed, and once errors in input data, i.e. surface data and prior information about basal properties are defined, the maximum a posteriori estimate and retrieval errors are uniquely determined.
We find that inversion of surface measurements of topography and velocity, combined with approximate estimates of mean ice thickness (h) and mean basal slipperiness, can be used to resolve spatial variations in basal topography down to about one ice thickness, with the exact number dependent on data errors and spacing. In comparison, the best obtainable spatial resolution of basal slipperiness variations is much more limited. For 1% data errors, for example, small amplitude spatial variations in basal slipperiness can only be resolved for wavelengths larger than about 50 mean ice thicknesses.
The number of resolved basal topography quantities per measurement site is close to unity down to a spacing between sites of oneh. For basal slipperiness the corresponding number is generally much smaller and only larger than 0.6 for spacing larger than tenh and for highly accurate measurements (data errors less than about 0.1% of mean values).
Accurate prior estimates of basal topography do not significantly improve the retrieval of variations in basal slipperiness. This does, however, not imply that measurements of basal topography are not useful in combination with surfaceto-bed inversion. Such independent information about basal topography can be expected to be useful for both model validation purposes and prior error covariance modelling.
The sensitivity of basal slipperiness retrieval to true basal topography increases sharply with decreasing prior errors in basal topography. Not allowing for some variation in basal topography when solving for spatial variations in basal slipperiness is equivalent to setting prior errors in basal topography to zero. In this case, estimates of basal slipperiness can be significantly affected by unmodelled errors in basal topography.