Estimating basal properties of glaciers from surface measurements : a non-linear Bayesian inversion approach

Estimating basal properties of glaciers from surface measurements: a non-linear Bayesian inversion approach M. J. Raymond and G. H. Gudmundsson Versuchsanstalt für Wasserbau, Hydrologie und Glaziologie (VAW), ETH Zürich, 8092 Zürich, Switzerland British Antarctic Survey, High Cross, Madingley Rd., Cambridge CB3 0ET, UK Received: 17 November 2008 – Accepted: 9 January 2009 – Published: 5 February 2009 Correspondence to: G. H. Gudmundsson (ghg@bas.ac.uk) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
The goal of geophysical inverse problems is to make quantitative inferences about Earth characteristics from indirect observations (e.g., Gouveia and Scales, 1998).Estimating basal properties of glaciers from surface measurements is an example of such an inverse problem.In this paper, we use a probabilistic Bayesian inference approach (e.g., Rodgers, 2000;Tarantola, 2005) to estimate bedrock topography and basal slipperiness from surface velocities and surface geometry.In Bayesian inference, a priori information about the basal properties is expressed as a probability density function and combined with the surface measurements to give a posteriori probability distribution describing the final uncertainty of the estimate.The solution of the inverse problem, i.e. the a posteriori probability distribution, provides an ensemble of solutions from which we single out the most likely one corresponding to the maximum of the a posteriori probability (MAP estimate).Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
The forward function describing the relationship between basal conditions (bedrock topography and basal slipperiness), and the observations (surface topography, surface velocities, rates of elevation change), is solved numerically with a non-linear planestrain finite-element model.A posteriori probability distribution for the system state is optimized via a nonlinear Gauss-Newton procedure to find the maximum a posteriori probability (MAP).
A key issue in the derivation of the solution involves determining the sensitivity of surface fields to perturbation in basal quantities.The Fr échet derivatives of the forward model can, in principle, be evaluated numerically.However, as the computational times involved in doing so are typically long in comparison to one forward model computation, it is, whenever possible, preferable to evaluate the derivatives of the forward function algebraically (Rodgers, 2000).Here we approximate the forward model derivatives using analytical transfer functions (Gudmundsson, 2003).These transfer functions describe the effects of small-amplitude perturbations in basal properties (bedrock profile and slipperiness) on surface fields in the case of Newtonian rheology and linear sliding law.
The numerical forward model solves for non-linear rheology, non-linear sliding law, and finite-amplitude basal perturbations.The transfer functions are, thus, only approximations to the actual forward model Fr échet derivatives.It is far from clear that using the analytical transfer functions in this context will result in a usable inverse model.What is clear, however, is that if this does work, the resulting improvements in computational efficiency are large.The main focuses of this paper is to determine whether this approximation is adequate in situation commonly encountered in glaciology.We do so by systematically constructing synthetic data sets were, to a varying degree, the assumptions of the analytical theory are broken.Thus, we start by using linear rheology and moderately strong basal perturbations, and then in a step-wise fashion introduce further sources of non-linear behavior.
We focus on situations similar to those on active ice streams.The method proposed here differs in a number of ways from previous inversion methods developed and applied to ice streams.Thorsteinsson et al. (2003) used the analytical small-amplitude Introduction

Conclusions References
Tables Figures

Back Close
Full solutions by Gudmundsson (2003) as a forward model for least-squares inversion of data from Ice Stream E, West Antarctica.In doing so the assumption was made that non-linear effects were not strong.A novel aspect of the method used by Thorsteinsson et al. (2003) was the simultaneous inversion for both basal topography and basal slipperiness.Gudmundsson (2006) suggested using a formal Bayesian inversion method instead of the least-square approach used by Thorsteinsson et al. (2003).MacAyeal (1992) and MacAyeal et al. (1995) applied control theory to determine the basal shear stress under ice streams using surface velocity data, ice thickness and surface elevation.Joughin et al. (2004) modified the method to yield a direct inversion for the basal stress corresponding to a weak plastic bed.These inversion procedures use forward models that solve reduced set of the Stokes equations.Another interesting approach to surface-to-bed inversion can be found in Truffer (2004) who inverted a linearized one-dimensional forward model to calculate the basal velocity of a valley glacier.Further examples of inversion of surface observations to determining basal conditions under glaciers can be found in e.g.Van der Veen and Whillans (1989) and Vieli and Payne (2003).
The purpose of this paper is to introduce and test the suitability of a non-linear Bayesian inference approach to determine the bedrock undulations and basal lubrication under ice streams from observations of surface topography, and horizontal and vertical surface velocity.An identical Bayesian inference approach, but for linear media and small basal amplitudes, has been presented in Gudmundsson and Raymond (2008).
The inverse procedure is applied along flow lines were transverse effects can be ignored.We use noise-degraded synthetical surface data generated with a forward numerical model.Finite-amplitude effects on rates of convergence are first investigated for linear rheology, and subsequently for non-linear rheology.The influence of nonuniform englacial temperature and non-linear sliding law on retrieval are also examined.The case studies presented allow us to explore the performance of the proposed inverse procedure and in particular to assess the practicality of approximating the Fr échet Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc

Printer-friendly Version
Interactive Discussion derivatives of the forward function using analytical small-amplitude solutions.We determine to what extent the inverse procedure converges to the true solution and how many iterations are needed.The structure of the paper is as follows.We start by describing in Sect. 3 the numerical forward model.We then introduce the non-linear Bayesian inference method in Sect. 4. The results of the proposed inversion method are presented in Sect.5, where the suitability of the method is discussed.

Notation
Vectors will be denoted by bold face italic letters (e.g.d) and matrices by bold uppercase letters (e.g.C).Surface measurements are available at discrete points, and we denote the set of all available surface quantities as the measurement vector d.The measurement vector d = (s, u, w ) T consists of surface topography s, horizontal velocity u and vertical velocity w .The basal properties to be estimated are assembled into one system state vector m.The vector m = (b, c) T contains the basal topography b and basal slipperiness c.The superscript T means transposition, here to column vectors.The subscript prior denotes a prior estimate, while a hat (e.g.m) indicates a maximum a posteriori estimate (MAP).

Forward model
The relation between basal properties and surface data can be written as (1) We refer to the function g as the forward model but we also use the term "forward function" when referring to g. the sliding law is non-linear, and because the surface reacts in a nonlinear fashion to finite amplitude basal perturbations.We use as forward function a numerical flow model model that allows us to deal with all types of nonlinearities mentioned above.The numerical model is two-dimensional finite element model that uses four-node isoparametric and quadrilateral Hermann elements.A mixed Lagrangian-Eulerian approach is employed in determining the position of the steady-state surface (Leysinger Vieli and Gudmundsson, 2004).The numerical model solves the full equilibrium equations (i.e.Stokes equations with acceleration terms set to zero), together with the mass-conservation equation for incompressible ice.These equations read σ i j,j = − ρg i and v i ,i =0, respectively, where σ i j are the components of the Cauchy-stress tensor, ρ is the ice density, g the acceleration due to gravity and v i are the components of the velocity vector v = (u, w).The glacier geometry corresponds to a uniformly inclined plane slab of constant thickness on which perturbations in bed and surface topography are superimposed.Figure 1 illustrates the problem geometry.
The coordinates are (x, z), where x is parallel and z perpendicular to the mean slope.The equation z=s(x, t) defines the surface and z=b(x) the base of the glacier.The constitutive law is Glen's flow law, extended, following Hutter (1983), with a linear term to avoid the singularity in viscosity as the deviatoric stress goes to zero (2) In this equation, A is the rate factor, n the stress exponent, ˙ i j , σ i j and τ are the strain rate, the deviatoric stress tensors and the effective shear stress, respectively.The parameter τ 0 is the crossover stress at which the linear and power terms contribute equally to the total strain rate.This parameter has been introduced only in model runs where the model geometry did not give rise to a sufficiently big mean longitudinal strain rate to avoid a large effective viscosity at the surface.Values for the rate factor for temperate ice are taken from Paterson (1994), and the dependency on temperature follows Smith and Morland (1981) Boundary conditions along the bed are specified by a sliding relation of the form where u b is the sliding velocity tangential to the bed, c(x) the sliding coefficient, and τ b the bed parallel shear stress.Basal sliding is introduced in the finite-element model by adding a uniform thin layer of different viscosity to the base of the glacier such that relation Eq. 3 is fulfilled.
The ice surface (z=s(x, t)) is stress-free and evolves with time according to the kinematic boundary condition until steady-state is reached.The kinematic boundary condition reads where s(x, t) describes the surface elevation, t is the time, and u and w are the horizontal and vertical velocity components, respectively.In model runs accumulation and ablation are not taken into account, but doing so is a straightforward modification.The kinematic boundary condition (Eq.4) is integrated forward in time with an semi-implicit Crank-Nicholson scheme.To speed up the evolution of the free surface toward steadystate, we initialize the computations with the analytical steady-state surface profile for linear rheology and respective basal perturbations (Gudmundsson, 2003).Periodic displacement boundary conditions are imposed along the upstream and downstream glacier model boundaries.
The size of the mesh in the x-direction follows a Gaussian repartition centered around the middle of the prescribed perturbations where we expect the largest deformations.This allows to reduce the total amount of elements and computational time.The results of the numerical model have been validated by comparison with relevant perturbation theories (Raymond and Gudmundsson, 2005).Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion

Perturbed fields
We define perturbations in boundary data and all field variables as the difference between the value of the variable in question at some given point and its spatial averaged mean value.For the basal topography, for example, we write where b(x) is the bedrock profile, b the mean bedrock elevation, and ∆b(x)=b(x) − b the bedrock perturbation.Similarly, the function describing the basal slipperiness c(x) is written as where c is the mean basal slipperiness, and c∆c(x) is the basal slipperiness perturbation.Perturbations in bedrock topography are referred to in the following as b perturbations, and perturbations in basal slipperiness as c perturbations.
Similarly to the basal perturbations, the steady-state surface topography (s), surface horizontal velocity (u), and surface vertical velocity (w) are partitioned as s(x)= s+∆s(x), u(x)= ū+∆u(x), w(x)= w+∆w(x), respectively.Collectively ∆u, ∆w and ∆s are the surface perturbations and ∆b and ∆c the basal perturbations.For reasons of notational compactness we will sometimes refer to s, u and w as the surface perturbations and to b and c as the basal perturbations.

Non-linear Bayesian Inversion
We perform a Bayesian inverse calculation to determine both bedrock topography and basal slipperiness from surface topography and surface velocities data.In Bayesian inference, the notion of knowledge and uncertainty about data and system state is expressed in terms of probability density functions (pdf's).The solution of the inverse problem is a posteriori probability distribution P (m|(d, m prior )) for the system state m, The denominator of Eq. ( 7) is independent of the system state and dose not affect the position of the maximum of the conditional probability on the left hand side of Eq. ( 7).
Hence, the a posteriori distribution reduces to the product of two terms, i.e. the a priori distribution P (m|m prior ) for the system state m and the likelihood function P (d|m).The likelihood function measures the probability of observing the data d if the system state was m, while the prior distribution incorporates prior information that is known independently of the measurements.As an example, the bed topography could be known independently from radio-echo sounding measurements.The a priori information may also arise from theoretical considerations (e.g., bedrock perturbation must be smaller than ice thickness), or some expectations (e.g., basal slipperiness not negative).Equation ( 7) is general.In this study, we assume that both data and modeling uncertainties can be described by Gaussian distributions.The a priori probability density function is therefore on the form where C M is the a priori covariance matrix describing the uncertainties in the prior system state.The likelihood function is given by

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
Here, C D is the covariance matrix of the noise in the data and g(m) is the forward modeling operator encapsulating the relevant physics in the relation between surface data d and system state m as described in more detail in Sect.3. Defining the cost function by and substituting Eqs. ( 8) and ( 9) into Eq.( 7) we obtain We solve Eq. ( 10) for the minimum of J(m) corresponding to the maximum of the posterior probability distribution P (m|(d , m prior )), that is, we single out the most likely system state m from the ensemble described by the pdf.This solution is referred to as the maximum a posteriori solution (MAP).Because the forward modeling operator g(m) is non-linear, there is no explicit solution to Eq. ( 10).We therefore perform a nonlinear optimization to find the maximum a posteriori solution m that maximizes P (m|(d, m prior )).To find the minimizers of J(m) we equate the derivative of Eq. (10) to zero. Defining where K(m)=∇ m g(m) is the Fr échet derivative matrix, the solution of the optimization problem is given by φ( m)=0.
The value of m is found using Newton's method via the iteration where the subscript i denotes the i -th iteration, the inverse is a matrix inverse and ∇ m φ(m) is the second derivative of the cost function, also called its Hessian.Equation ( 13) involves both the first derivative K(m) and the second derivative ∇ m K(m) of the forward model, whose resulting product is small in the moderately non-linear case (Rodgers, 2000).Ignoring this term and substituting Eqs. ( 11) and ( 13) into the Newtonian iteration Eq. ( 12) gives the m i +1 iteration according to the Gauss-Newton method, namely where The Fr échet derivatives of the forward model K i , are approximated using linear transfer functions, i.e.
The transfer functions T are analytical solutions for linear rheology describing the effects of small-amplitude variations in bed topography and basal slipperiness on surface fields (Gudmundsson, 2003).The transfer functions T have a two letter suffix.The first suffix denotes the effect and the second one the cause.T sb describes a change in surface topography caused by a perturbation in bedrock topography, whereas T uc describes a change in surface-parallel velocity caused by a spatial variation in basal slipperiness.Figure 2 shows examples of analytical transfer functions as functions of the wavelength for both b and c perturbations.When using the transfer function formulation for K i , the inversion is most easily done in Fourier space.In Eq. ( 14 transformed and H the Hermitian transpose.The transpose T in Eq. ( 14) is substituted with the Hermitian transpose H.

Data uncertainties
The covariance matrix for the noise in the data C D is defined as The matrix C D is a block diagonal matrix consisting of the matrices describing the uncertainties in the surface topography C s , horizontal velocity C u and vertical velocity C w along the main diagonal.The off-diagonal blocks are zero matrices, since no cross-correlation between errors in surface topography, horizontal and vertical velocity is considered.The covariance matrix for the noise in surface topography C s takes the form where N is the number of discretization points.The matrix element σ 2 s i s j is the covariance of s i and s j , with 1≤i , j ≤N.If the noise in the data is uncorrelated, the corresponding covariance matrix is of diagonal form.C u and C w are of the same form but with the indices s replaced by u and w, respectively.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
The covariance matrix for m prior is defined as where C b and C c have the same structure as C s .No a priori cross-correlation between the a priori estimates of basal topography and basal slipperiness is considered.
The correlation between off-diagonal elements is assumed to follow a Gaussian distribution and the elements of C b are therefore given by σ b describes the variance of the fluctuation about the mean m prior of the Gaussian probability density that characterizes the uncertainty of the a priori estimate, l b is the correlation length and x the position.A corresponding expression is used for basal slipperiness a priori covariance matrix (C c ) with l c denoting the correlation length.

First guess
To start the Newtonian iteration we need an initial guess for m i =0 in Eq. ( 14).The initial guess generally corresponds to the a priori values for the model parameters m prior .As stated above, we consider the case where m prior = 0. To define m i =0 , we assume that the relationship between basal and surface properties is linear (i.e. the ice rheology is linear and the amplitude of the basal perturbation is small) and can be completely described using the perturbation theory of Gudmundsson (2003).Thus, the forward relation Eq. ( 1) simplifies to d =Km + , where K is a matrix of transfer functions.As all pdf's are still Gaussian, the cost function is of same form as Eq. ( 10) but with g(m) replaced by K.
Hence, we start the optimization of the objective function J(m) by setting m i =0 = mlinear .

Convergence
The iteration is continued until either the convergence test is fulfilled, or maximum number of iteration has been exceeded.The above criterion is based on the fact that 3×N corresponds theoretically to the expected value of J( m) (Tarantola, 2005).iv Inversion step: Determine incremental corrections to the bedrock profile and the distribution of the basal slipperiness using ( 14).Return to step (ii) Note that the analytical transfer functions are used not only in the initialization step (i) but also in step (iv).

Model experiments
The main focus of this paper is to evaluate the performance of the inverse method described above in situations typically encountered on active ice streams.The key novel aspect of the method is the use of the analytical transfer functions in approximating the forward model Fr échet derivatives.This approximation renders the method tractable on current generation of computers.On the other hand, because the transfer functions are based on small-amplitude perturbations theory, this aspect of the method is also the one which might most severely restrict its applicability.We start by inverting for basal properties for cases where we expect the method to work, i.e. for linear rheology, linear sliding law, and small perturbation amplitudes.Subsequently we introduce large perturbation amplitudes, non-linear ice rheology, and non-linear sliding.Only a subset of the experiments performed are presented here.A more detailed description can be found in Raymond (2007).

Synthetic data
All synthetic surface data sets (i.e.surface topography, horizontal and vertical velocity) were generated with the finite-element forward model described in Sect.The synthetic surface data set was corrupted by uncorrelated Gaussian noise.The forward relation Eq. ( 1) hence modifies to where are uncorrelated measurement errors.
The surface data was assumed to be available at equally spaced location along the profile.The surface data was interpolated to the nodes of the non-equidistant finite-element mesh.This interpolation introduces errors that, in general, are spatially correlated, and the covariance matrix describing the surface data errors is therefore no longer of diagonal form.Data was interpolated using the Best Linear Unbiased Estimate (BLUE) (e.g.Kitanidis, 1997).The method makes use of the experimental variogram to describe the spatial variability of the measurements.The variogram was approximated with an isotropic Gaussian model.Key aspect of this interpolation method is that it delivers full covariance matrices C s , C u , and C w , describing the combined effects of noise in the original surface data and the errors and correlations introduced by the interpolation.

Dimensions
All results are presented in non-dimensional form.The dimensional variables entering the problem are substituted through scalings by non-dimensional variables (Raymond and Gudmundsson, 2005).The spatial variables (x and z) are scaled with the mean ice thickness h and the velocity components (u and w) with the mean surface-parallel deformational velocity ūd .The stresses and the pressure are scaled by the mean basal Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion shear stress τb .For a given mean surface slope, mean ice thickness, and stress-law exponent, the rate factor is scaled so that the mean forward deformational velocity is equal to unity.It follows that in non-dimensional form, the rate factor A equals to A = (n+1)/2.When inserted into the sliding law these scalings give a non-dimensional basal slipperiness C(x)=c(x)τ (0)m / ūd where as before c(x) is the dimensional slipperiness and C(x) is its non-dimensional counterpart.It follows that C = ūb / ūd , where C is the mean non-dimensional slipperiness (C(x)= C(1 + ∆C)) and ūb and ūd are the mean basal and deformational velocities, respectively.C is therefore equal to the slip ratio which is defined as the ratio ūb / ūd .
5.3 Linear rheology, linear sliding law, and small-amplitude perturbations Figure 3 illustrates a synthetic surface data set generated with the numerical forward model.The data represents surface response to a combined perturbation in both bedrock profile and basal slipperiness.The amplitudes of both basal disturbances are small, or 5% of their respective mean values, and the perturbations are Gaussian shaped with standard deviation of 12 mean ice thicknesses.Mean surface slope was 0.5 degrees and mean slip ratio equal to 200.Further information on basal parameters is given in Table 1 and in the figure caption.
The surface data were corrupted by uncorrelated Gaussian noise with zero mean and standard deviation equal to 5×10 −4 for the surface topography, 10 −3 and 3×10 −3 for the horizontal and vertical velocity components, respectively.This level of noise is realistic for situations where surface measurements are done using GPS technology.Surface measurements were equally spaced with an interval equal to two ice thicknesses.The interpolation of surface data one the surface nodal points of the finiteelement mesh was done with the method of best linear unbiased estimator (BLUE).The surface data covariance matrix used in the inversion is the corresponding BLUE estimate of interpolation errors.
The a priori estimate was no perturbation in either basal topography or basal slipperiness with large unknowns.The a priori covariance matrix was a Gaussian with 197 Introduction

Conclusions References
Tables Figures

Back Close
Full diagonal elements corresponding to a 25 % error in the respective mean values of basal slipperiness and ice thickness.The a priori correlations lengths for both the c and the b perturbation were equal to 10 mean ice thicknesses.
Figure 4 shows the inferred bed topography and basal slipperiness distribution using the data shown in Fig. 3.The first guess, the subsequent iteration, and the maximum a posteriori solution (MAP) obtained after two iterations are all shown.For comparison the true basal disturbances are depicted using solid lines with circles.As can be seen in Fig. 4, the first guess (blue lines), obtained by assuming a linear relationship between basal and surface properties, already resolves both basal perturbations quite accurately.This is to be expected since the ice rheology is linear and perturbation amplitudes fairly small (5%).Nevertheless, the subsequent non-linear optimization step significantly improves on the initial estimate, and the final maximum a posteriori solution ( m) obtained after only two iterations is almost identical to the true value of the basal disturbances (m).
Figure 5 shows the residuals between observations and FE-model predictions for, (a) the surface topography ∆s, (b) the horizontal velocity ∆u and, (c) the vertical velocity ∆w for iteration number 0 and 2. The residuals are defined as ∆d i =d − g(m i ) where ∆d i is the vector of residuals and g(m i ) the forward finite-element model prediction for the system state m i at iteration numbers i =0, 2. The optimization procedure was continued until the forward model results agreed to observed data to within one standard deviation (dotted lines).As can be seen, only the horizontal velocity is poorly predicted by the FE-model by the first guess.This is almost definitely because the slipperiness perturbation, which has a much larger effect on horizontal velocity component than the vertical one and negligible effect on surface topography, is not resolved in full detail in this initial step.Note that even in the final MAP estimate the b perturbation is better resolved than the c perturbation.As we will see below, this is generally the case.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion 5.4 Linear rheology, linear sliding law, and large amplitude Gaussian perturbations Figure 6 shows the inferred bedrock and basal slipperiness perturbations for the Gaussian peak distributions as used in the example above, but now for perturbation amplitudes of 20% instead of 5%.All other parameters are the same as in previous 5% amplitude example.
Whereas for 5% amplitude perturbation it would have been surprising had the method not converged in a few iterations, we now have basal amplitudes that can be expected to be too large to be well retrieved in the first initialization step (Step (i)).In addition, the forward model Fr échet derivatives used in step (iv) are now considerably less accurate.As expected on the basis of this, the first guess retrieval is poorer than it was in the 5% amplitude case.The initialization step creates a peak in the basal slipperiness at the right position (x= − 10).However, it also produces a negative peak at x=10 with a considerable amplitude of more than −0.1.
Encouragingly, a positive perturbation of smaller amplitude appears at the same location in next iteration step (Fig. 7).In the following iterations the retrieval converges quickly to the true solution.Only two more iterations are needed as compared to the 5% amplitude case, and only a total of four iterations is needed to compute the the MAP estimate.

Non-linear ice rheology, linear sliding law, and small-amplitude Gaussian perturbations
We now turn our attention to the case of non-linear ice rheology, keeping for the time being the sliding law linear.Calculations were done for a stress exponent of n=3 (see Eq. 2).The 5% amplitude experiment with the two Gaussian shaped bedrock and slipperiness perturbations described above was repeated, this time using n=3 rather than n=1.Despite the introduction of non-linear ice rheology the inversion algorithm not only converged but did so in only 7 iterations.Further details of this experiment and other Introduction

Conclusions References
Tables Figures

Back Close
Full small-amplitude experiments using bump-shaped sigmoidal functions can be found in Raymond (2007).

Non-linear ice rheology, linear sliding law, and large-amplitude Gaussian perturbations
We now consider the inversion of noise-degraded synthetic data generated for nonlinear rheology and relative large amplitudes of the basal perturbations of 20%.Both rheological and finite-amplitude nonlinearities are thus present and affect the surface response.The cross-over stress was τ 0 =0.3 τb .
The widths of the basal perturbations, surface data errors, slip ratio, mean ice thickness, and a priori estimates are all the same as before and listed in Table 1.However the mean surface slope is now α=0.2 • instead of α=0.5 • .The results are not dependent on the exact value of the surface slope.
The maximum a posteriori solution, computed using 11 iterations, is shown in Fig. 8 along with the true basal disturbances (black lines with circles).As can be seen, the iteration procedure converged to the correct solution and the MAP reproduces quite accurately the prescribed basal disturbances.The initialization step is somewhat inaccurate (see Fig. 8) and the amplitudes of the initial guess for the b and the c perturbations had to be forced to remain smaller than 30% of the respective mean values.In comparison to the linear-rheology case twice as many iterations are needed for convergence.Nevertheless, this is a very good result as it shows that approximating the forward model Fr échet derivatives with the analytical transfer function in the optimization procedure (Eq.14) is possible in this context even when both nonlinear rheology and finite-amplitude effects are combined.
Looking at the residual (Fig. 9) shows that the iteration is in a the form of a damped oscillation toward the final solution.This is seen in both the b and the c perturbation fields.After iteration 3, the oscillating behavior disappears.Introduction

Conclusions References
Tables Figures

Back Close
Full In all previous experiments the sliding law was linear (i.e.m=1).In the following, we will consider the inversion of noise-degraded synthetic data generated for non-linear ice rheology and non-linear sliding law, i.e. m=3 in Eq. 3. The surface data were computed for Gaussian peak distributions with perturbation amplitudes of 5%.In this experiment, the surface slope α=0.2, the slip ratio C = 200, τ 0 =10 kPa, the mean ice thickness h=1000 m and the deformational surface velocity u d =4.03 m a −1 .The inversion procedure converged to the correct solution in 10 iterations.Further details of this experiment can be found in Raymond (2007).
5.8 Non-linear rheology, non-linear sliding law, and large-amplitude Gaussian perturbations The following experiment considers amplitudes of the basal perturbations of 20%.All others parameters are the same as in the previous 5% amplitude non-linear sliding law case.Figure 10 shows the inferred bedrock and basal slipperiness distributions (red lines).As for non-linear ice rheology, linear sliding law and large-amplitude perturbations (Sect.5.5, see also Fig. 8), the amplitudes of the initial guess for the basal slipperiness had to be forced to remain smaller than 50% of the respective mean values (blue lines in Fig. 10) by imposing a corresponding constrain on the a prior estimate.Interestingly, the form of the initial basal slipperiness distribution is similar to the corresponding linear-sliding case (Sect.5.5), but with a larger amplitude.The MAP solution is obtained in 14 iterations and reproduces quite accurately the prescribed perturbations.As compared to the corresponding linear-sliding case, only 3 more iterations are needed for convergence.Figure 11 shows the residuals between observations and FE-model predictions.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion 5.9 Temperature-dependent non-linear ice rheology We now introduce two modifications to the experimental setup, (1) the rate factor is no more constant across the thickness, and (2) we use basal perturbations corresponding to bump-shaped sigmoidal functions instead of the Gaussian used above.The vertical variation in rate factor corresponds to a linear variation in ice temperature from −25 • at the surface to 0 • at the base.The rate factor A is then determined using the double exponential fit between A and temperature derived by Smith and Morland (1981).The surface-parallel deformational velocity amounts in this case to u d =1.27 m a −1 .
The basal disturbances are shown in Fig. 12.The parameters are: mean ice thickness h=300 m, surface slope α=1 • , slip ratio C = 40, cross over stress τ 0 = 0 kPa, and 5 % bump-shaped perturbations in b and c (see also Table 2).Surface data errors and a priori estimates are the same as in all previous experiments.
As can be seen in Fig. 12 the (linear) initialization step (blue lines) already gives a good approximation to the true bump-shaped perturbations in b and c.In comparison to the retrieval of the b perturbation, the c initial estimate is, however, less accurate and the rate of convergence slower (see Fig. 12b) and a total of 14 iterations is needed for convergence.In contrast to Fig. 8, the estimates converge toward the solution without oscillations and the residuals (Fig. 13) diminish continuously toward the noise-level.
In Raymond (2007) a corresponding experiment using isothermal conditions is described.Comparisons reveals that about four times as many iterations are needed than for the isothermal case to fit the observed data down to the noise level.The final MAP model estimates are, however, equally good showing that only the rate of convergence but not the quality of the solution suffers when the rate factor varies across the depth.

Forward model parameter error
The output of the forward model depends not only on the basal boundary conditions, that is on the form of the c and the b perturbations, but also on a number of forward model parameters.There is a total of six such forward model parameters: the stress Introduction

Conclusions References
Tables Figures

Back Close
Full exponents n and m, mean ice thickness h, and mean surface slope α, the mean slipperiness C, and the rate factor A. In all experiments described above it has been assumed that the values of these model parameters are known with perfect accuracy, whereas in a real situation these parameters are always known within some error bounds.
The general approach to this problem would be to include the forward model parameters as elements of the system state and formally invert for these parameters in a manner identical to the inversion for b and c.However, a number of these parameters are either rather well constrained through laboratory measurements (A and n) or can be expected to be rather easily determined from the data (mean surface slope, mean ice thickness).It seems likely that the mean basal slipperiness ( C) and the stress exponent of the basal sliding law (m) are the forward model parameters most likely to be estimated incorrectly.We do not attempt here to give a general description of errors in retrieval due to forward model parameter errors.We limit the discussion here to giving an example of a retrieval where the mean value of the slipperiness is only known approximately.
Let us consider the synthetic surface data generated for bump-shaped perturbations with linear temperature profile increasing from −25 • at the surface to 0 • at the bottom (see the case for temperature-dependent ice rheology in Sect.5.9).The true deformational velocity amounts to 1.27 m a −1 and the corresponding slip ratio is C=40.As the temperature profile is not known, we start by estimating a linear temperature profile with surface temperature −20 • and bottom temperature 0 • .The mean surface deformational velocity u d is estimated to be 1.45 m a −1 using the standard temperature-dependent flow law for ice.The slip ratio is then estimated from the mean longitudinal surface speed to C=35.These estimated values of the forward model characteristics are subsequently introduced in the inverse calculations as well as in the forward finite-element model.
Figure 14 shows a comparison of the maximum a posteriori model as obtained with the true model characteristics and with the modeling errors.The true basal properties are also shown for comparison.In both cases, the basal topography perturbation Introduction

Conclusions References
Tables Figures

Back Close
Full is completely recovered.The basal slipperiness perturbation is, on the other hand, resolved less accurately for the case with the wrong estimated model characteristics than with the true ones.The inversion procedure has been stopped after 7 iterations for the wrong model and after 14 iterations for the true one.The value of the maximum a posteriori solution J( m) at the minimum for the case with the wrong model characteristics is 8 times larger than for the true model and also about 8 times bigger than the theoretically expected value of J( m)=3×N at the minimum.Thus, if we were to select a forward model characteristics set best describing the observations, we would choose the one for which the value of the maximum a posteriori solution J( m) is closest to the theoretical one.The MAP solution for the wrong model characteristics would obviously be discarded.

Conclusions
We have shown how a non-linear Bayesian inference approach can be used to simultaneously determine both basal topography and basal slipperiness from surface measurements of velocity and topography.We empathize the fact that we do not produce one single solution to the inverse problem considered.Rather, we determine the maximum a posteriori solution together with the a posteriori error covariance matrix.
The availability of such an error estimate greatly facilities any quantitative analysis of the results.Furthermore, the method does not require any prior smoothing of input data.
We have shown that the Fr échet derivatives of the forward model can be adequately approximated by small-amplitude analytical solutions (Gudmundsson, 2003) for the method to converge.Our key result is that this remains the case even when the problem is strongly non-linear.This result is of considerable practical value as this approximation greatly enhances the numerical efficiency of the method by sparing the time-consuming numerical evaluation of model derivatives.In fact, we find it difficult to envision that an inverse method employing the full Stokes equations in the forward step Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion could be of any practical value unless a good way of estimating the Fr échet derivatives is found.The proposed method is capable of dealing with both non-linear finite-amplitude effects and rheological nonlinearities.In all case studied, the inversion procedure converged quickly and to the correct solution and in only a few iterations, with the exact number of iterations needed being dependent on the type and magnitude of non-linear effects.Introduction

Conclusions References
Tables Figures

Back Close
Full   Full ) all vector components, i.e., surface fields, a priori and basal perturbations are therefore transformed to frequency space by the relation Fs where F is the unitary discrete Fourier transform matrix and s the vector to be transformed.The covariance matrices for the data and model parameters C D and C M are transformed to the Fourier space by the relation FCF H where C is the matrix to be Introduction

4. 4
Inverse procedureIn summary, the different steps involved in the iterative optimization by which the objective function J(m) is minimized, are: i Initialization step: Define a first guess for m i =0 assuming everything can be described perfectly by the analytical transfer functions valid for small-amplitudes and Newtonian rheology.iiForward step: Calculate the steady-state surface response g(m i ) for the given bedrock and the distribution of the basal slipperiness with the non-linear forward finite-element modeliii Convergence test: Test for convergence using Eq.(20).Once the stopping criterion is satisfied, stop the iteration procedure, else 3. The perturbations in bedrock and basal slipperiness correspond either to Gaussian peak distributions ∆b(x) = a b exp − (x − x b ) 2 expressions for the c perturbations.Here a b is the basal amplitude, x b the center of the perturbation and the standard deviation σ b defines the width of the basal perturbation.
-linear rheology, non-linear sliding law, and small-amplitude Gaussian perturbations

Fig. 1 .
Fig. 1.Illustration of the problem geometry and coordinate system.Gaussian-shape bedrock perturbation and corresponding surface reaction are shown as a black line.The dashed lines show the undisturbed glacier geometry.

Table 1 .
Values of the model parameters for all Gaussian perturbation experiments.Here n is the stress exponent in Glen's flow law, a b and a c are the amplitudes of the bedrock and basal slipperiness perturbations, respectively, α is the mean surface slope and H the mean ice thickness.In all Gaussian experiments, the center of the bedrock and basal slipperiness perturbations are located at x b =10 and x c = − 10, respectively.The widths of the basal perturbations are described by the standard deviation σ b =σ c =12 mean ice thicknesses.The mean basal slipperiness is C = 200.

Table 2 .
Values of the model parameters for the bump-shaped sigmoidal perturbation experiments.Here a b and a c are the amplitudes of the bedrock and basal slipperiness perturbations, respectively, α is the mean surface slope and H the mean ice thickness.The stress exponent in Glen's flow n=3 in both experiments and the mean basal slipperiness C=40.