2.1. Kinematic equations

The approach presented here is based on a kinematic consideration of the firn volume; the equations for conservation of energy and momentum are therefore not taken into account. In general, the distribution of any tracer in a medium can be described by an advection–diffusion equation. Details of the tracer transport and formulation in ice sheets are discussed extensively by Reference Clarke, Lhomme and MarshallClarke and others (2005). In our case, the corresponding tracer is depositional age A = A(r, t), a non-diffusive property, which obeys

(1)

where v = (u, w) = v(r, t) is velocity, ∂ _t denotes partial derivative with respect to the subscript variable (here time t) and all calculations are carried out in two-dimensional (2-D) space r = (x, z) (z positive and increasing downward). See the Appendix for conventions and a list of notation.

Equation (1) is sometimes referred to as the age equation (e.g. Reference Hindmarsh, Leysinger Vieli, Raymond and GudmundssonHindmarsh and others, 2006). The righthand side represents a source term, which is responsible for the actual aging of the firn with time.

The second governing equation is the conservation of mass,

(2)

where ρ = ρ(r, t) is density. These two equations form the fundamental system of linear equations used in the forward problem.

2.2. Assumptions and boundary conditions

A number of assumptions are made for the sake of simplicity; however, they do not reduce the general applicability of the inverse-problem formulation. The considered firn volume extends from the surface z = 0 to an arbitrary depth z = z _max. The density distribution is assumed to be laterally homogeneous and independent of time, i.e. ∂ _xρ = ∂ _yρ = ∂_t ρ = 0 (Sorge’s law), but depth dependence is maintained (∂ _zρ ≠ 0). This assumption is well justified on a regional scale for ice-sheet plateaus (e.g. Reference FrezzottiFrezzotti and others, 2004; Reference Richardson-NäslundRichardson-Näslund, 2004; Reference Rotschky, Eisen, Wilhelms, Nixdorf and OerterRotschky and others, 2004; Ar-cone and others, 2005), but has to be considered with care on cold alpine glaciers.

Note that the depth dependency of density is a prominent deviation from the incompressibility assumption often used in ice-sheet modelling. Time dependence of Equations (1) and (2) is maintained in the prognostic forward model. The system of equations to be solved, however, is formulated in a time-independent way so that ∂_t (·) = 0 (where (·) denotes any term to be differentiated) as the forward model produces a steady-state age distribution as output.

No forces appear in the above equations, simplifying matters such that the upper boundary can be taken as a horizontal surface, i.e. parallel to x. Position and direction of scalar and vector quantities therefore always refer to this surface. (Consider a radargram as an illustrative example. It contains records of the reflector depth with respect to the relative surface. A topographic correction is applied only during data processing.) The kinematic boundary condition at the surface is , where is the surface accumulation and ρ ₀ = ρ(z = 0) is the density at the surface. Additional constraints are introduced in section 4.2, primarily as prescribed velocity properties.

2.3. Prognostic forward model

The forward model runs under prescribed stationary allocations of density, horizontal velocity and accumulation on an ordinary grid, discretized with finite differences. It calculates the vertical velocity from the combined effect of accumulation at the surface, advection and densification, and yields the synthetic age–depth distribution. Starting from an initial laterally homogeneous, vertically increasing age distribution, the prognostic model runs in a transient mode until a steady state is reached, i.e. when the particles from the surface at t = 0 reach the edge of the domain. As a boundary condition, age is set to zero at the surface. At the inflow of the model domain, the horizontal age gradient is set to zero.

Details of grid parameters are listed in Table 1. The age–depth distribution constitutes the essential output, which is passed to the inverse problem. The prescribed horizontal velocities u ^ref and calculated vertical velocities w ^ref of the forward model are defined for all gridpoints. We later refer to them as the reference-velocity field (denoted by the superscript ‘ref’), against which the inferred velocity field (denoted by the superscript ‘est’) is compared.

Table 1. Simulation parameters

2.4. Linear system for inverse model

The time-independent forms of Equations (1) and (2) are

(3a)

(3b)

The discretization schemes for solving this linear system on a triplex-staggered grid (a grid consisting of three sub-grids shifted relative to each other) are adapted from Reference Fiadeiro and VeronisFiadeiro and Veronis (1982) and Reference WunschWunsch (1985). The input fields of age and density are prescribed on a rectangular grid, the A-grid, with a grid spacing of Δx and Δz in the x and z directions, respectively. The A-grid has I × K nodes. Corresponding indices for the gridded variables are i = 1, …, I for the horizontal coordinate (increasing downstream, left to right) and k = 1, …, K for the vertical coordinate (increasing downwards, top to bottom), as indicated in Figure 1a. The grid nodes representing u and w (u and w grid) are shifted by half the grid spacing in the horizontal and vertical direction, respectively, relative to the nodes on which the input parameters for age A and density ρ are prescribed (Fig. 1).

Fig. 1. (a) Unit-cell scheme of the numerical grid used for solving the linear system of Equation (3). (b) Scheme of the triplex-staggered numerical grid for I = K = 6. The uppermost row corresponds to the surface. Distance between nodes of similar types is Δx and Δz and between nodes of different types is Δx/2 and Δz/2 in the horizontal and vertical directions, respectively. The cross centred on the (⊗-node labeled A _2,4 represents the unit cell in (a) and strikes all nodes involved in the age equation for the A _2,4 node. Likewise, the cross labeled ρ _4,3 strikes all nodes involved in the conservation-of-mass equation for the ρ _4,3-node. Both equations can therefore only be solved for those A-nodes within the region bounded by the dashed line, referred to as solution domain. The (&-nodes on the corners are displayed for completeness, but not used in the inverse problem.

The application of staggered-grid differences to Equation (3) leads to a discrete system, which for a unit cell (Fig. 1a) can be expressed as

(4)

Detailed expressions of the staggered-grid differences and coefficients are given in the Appendix. As sketched in Figure 1b for the node labeled A _2,4, five A-nodes are involved in the discretized representation of the age equation for a single node. Consequently, the u_{i ,k} , w_{i ,k} for a unit cell always depend on the values of A and ρ at the neighbouring nodes. These values are contained in the c_{i ,k} -coefficients in Equation (4). The u_{i ,k} , w_{i ,k} therefore cannot be fully determined on the boundaries, but only within the dashed region shown in Figure 1. This region is termed the solution domain. This formulation has the advantage that no other specific conditions are necessary at the boundaries of the domain where the inverse problem is solved with SVD.

As can also be seen in Figure 1b, in each dimension, x and z, the total number n of nodes for unknown variables u and w differs. Within the solution domain, the number of variables u in a row (x direction) is . Analogously, the number of variables along a column (z direction) is . For the variable w, and . The total number of elements of each variable within the solution domain is and . Defining the vectors and matrix

(5)

allows us to set up the matrix equation

(6)

The variables p, q are merely indices of vector and matrix elements, to be distinguished from the coordinate indices i, k of the actual grid. The vector d represents the data in data space , and the vector v represents the model parameters in model space is the number of (known) equations and N is the number of unknowns, in our case the velocities within the solution domain. The relationship between model parameters and data is described by the model matrix M, sometimes referred to as the data kernel (Reference MenkeMenke, 1989, p. 9).

The reader might wonder how it is actually possible to define uncertainties of the data vector d, which contains only ones and zeros. For this particular inverse problem, the measurable quantities age and density appear on the lefthand side in the matrix elements of the data kernel. The uncertainty of the data vector is therefore a measure of how the uncertainties of the data kernel cause the vector on the righthand side of Equation (6) to differ from exactly ones and zeros, even for exact velocities v. This point is explored in more detail using a Monte Carlo-based approach in section 5.5.

3.1. Principles

The SVD of a matrix M is a generalization of the spectral decomposition of a square to a rectangular matrix. The spectral decomposition of a rectangular matrix always exists. Here we apply SVD to calculate the pseudo-inverse (or generalized inverse) of M, mainly following the notation of Reference WunschWunsch (1996). Any rectangular matrix M can be decomposed into a factorization of the form

(7)

where U and V are both unitary rectangular matrices, , and V ^T denotes the transpose of V. The generally non-square matrix contains the singular values (square root of eigenvalues) of M in decreasing order on the main diagonal, Λ _{p ,q} = δ_pqλ_p , with the Kronecker symbol δ_pq . The matrix V contains a set of orthogonal base-vectors of M, spanning the N-dimensional model (or solution) space, whereas the matrix U contains a set of orthogonal base-vectors spanning the M-dimensional data (or observation) space. The number R of non-zero singular values is the rank of M. If some singular values are zero or M ≠ N, one or more of the rows or columns of Λ must all be zeros. Those columns of U and V that are multiplied by zeros only can be dropped, thus reducing the matrices in Equation (7) to the expression

(8)

where the subscript R indicates the number of columns with and . is the square submatrix of Λ with non-vanishing singular values. It can be shown (e.g. Reference WunschWunsch, 1996) that is the pseudoinverse of M, which we use to solve Equation (6) for the unknown model vector. We obtain the solution

(9)

where is the inverse of Λ _R , i.e. with on the main diagonal (λ_p ≠ 0) and zeros elsewhere.

The above expressions for M, U and V define four spaces: the model range (column space of M); the model nullspace ; the data range (row space of M); and the data nullspace . Depending on the size of M, N and R, not all of these spaces need to exist (in the sense that they are not empty sets). Conditions for existence of these spaces, definition for over- and underdetermined systems of equations and combinations of these are listed in the Appendix.

If there is a data nullspace U ₀ (R < M), and if the data have components in it, then it will be impossible to fit the data exactly. This data mismatch between true data and estimated data, referred to as the residual norm, will then be different from zero. (As a norm we will use the L ₂ norm or Euclidean length of a vector, denoted by the operator ∥ · ∥; see Appendix for definition and further information.) On the other hand, if the model has components in the model null-space V ₀ (R < N) then it will be impossible to determine the model exactly: hence the term model nullspace. In that case, the model solution can be presented as a sum of the particular solution given by Equation (9), which contains only range vectors and solves Equation (6), and an arbitrary homogeneous solution V ₀ α which solves the homogeneous system of equations Mv = 0. The vector α contains (N−R) coefficients for the linear combination of the (N−R) column vectors of V ₀ in the model nullspace, about which the equations provide no information.

The SVD is related to the least-squares approach. All of the structure imposed by SVD is also present in least-squares solutions. One commonality is that the SVD simultaneously minimizes the residual and solution norms (e.g. Reference Scales, Smith and TreitelScales and others 2001, p.66). However, the SVD solution generalizes the least-squares solution to the case where the matrix inverses of M ^T M or MM ^T , the simplest forms, do not exist, for example if the system is not full rank (Reference WunschWunsch, 1996, 157f). An important advantage of the application of SVD and the interpretation of the solution is that only a single algebraic formulation is necessary for over-, under- or just-determined systems. The SVD provides its control over the solution norms, uncertainties and covariances through choice of the effective rank , which leads to the so-called truncated SVD, demonstrated in section 5.1. The truncated form makes a clear separation between range and nullspace in both solution and data spaces.

3.2. Resolution

A useful feature of the SVD is that it provides direct access to the resolution obtainable when mapping between model and data spaces (see Reference MenkeMenke, 1989, p. 62f; Reference WunschWunsch, 1996, p.165). The model resolution matrix, defined as

(10)

determines the relationship between the general solution and the particular solution. If no model nullspace exists (R = N), the general and particular solution are equal. Then T_V = I _N , the N × N-dimensional identity matrix, meaning that the model is completely resolved. If a nullspace exists, non-zero terms will appear off the main diagonal in Equation (10), so only averages of some model parameters can be resolved.

Analogously, the data resolution matrix

(11)

provides information on how well the observed data are estimated when the model solution obtained with the generalized inverse is used in the forward model to predict observable quantities.

Both resolution matrices are functions of the data kernel M, which contains the a priori information about the physical representation of the problem, i.e. by the time-independent Equations (3). When the problem is linear, resolution matrices depend on neither the model parameters v nor the data d.

3.3. Error covariance and uncertainty

Solving the inverse problem yields an estimate of model parameters, denoted v ^est, which are subject to uncertainties. Using the estimated v ^est in the forward problem, Equation (6) yields a prediction of the data vector d ^est which differs from the true data vector d by some residuals, denoted n = d−d ^est. The residuals can generally arise from two contributions: noise from errors in the measurement of data and inadequacy of the forward algorithm to describe the problem exactly.

The covariance C_vv of the estimated model parameters depends on the residual covariance R_nn (the second moment or covariance matrix of n; see Appendix for details). It can be shown to be (Reference WunschWunsch, 1996, p.143)

(12)

In the case of uncorrelated uniform variance of the data, Equation (12) simplifies to

(13)

The covariance of the model parameters arises from uncertainties present in the data and generates uncertainty in the coefficients of the model range vectors. Data covariance is thus mapped onto model covariance. To obtain the complete solution uncertainty P_vv of the model parameters, the influence of the missing nullspace contribution also has to be taken into account. According to Reference WunschWunsch (1996),

(14)

where R _αα is the second-moment matrix (or covariance matrix; see Appendix) of the coefficients _α of the model nullspace V ₀, forming the homogeneous solution V _{0 α} . The matrix R _αα may be entirely unknown, or an estimate from a priori information might be available. The uncertainty of the residuals P_nn follows from the variance of the estimated residuals about their mean (Reference WunschWunsch, 1996, p.117), which can be written as

(15)

The covariance Equation (12) of the estimated model parameters is very sensitive to small non-zero singular values. Solution variance can be reduced by choosing an effective rank to exclude small λ_p . Inspecting the singular-value spectrum of the data kernel enables one to choose an appropriate cut-off size for contributing singular values (Reference MenkeMenke, 1989, p.122). This artificial reduction of model- and data-space dimensions leads to rank deficiency and therefore lower resolution as well as increased dimensions of the nullspaces, but decreases model covariance. The choice of the effective rank therefore provides a means to trade-off variance and resolution, or solution norm and residual norm, respectively.

3.4. Scaling and weighting

Weighting is in general used to give more importance to certain observations than to others, mainly to correct for uncertainty. An undesired weighting effect occurs if a system consists of different physical equations, involving different physical quantities.

In our case, the conservation-of-mass and the age equation involve the quantities age and density. In the linear system Equation (6), the rows of M represent these equations. Their different physical origin leads to different norms of the row vectors (i.e. Euclidean length) of the matrix M.

To correct for this effect, we first perform row scaling of the matrix M by multiplying each row with the reciprocal of its row norm (see Appendix for details). This is carried out below by operations with the matrix W, which contains the row norms of M on its diagonal. Likewise, the column vectors of M have different norms. We therefore require column scaling after the row scaling is performed. This is done by operations with the matrix S. Performing row scaling first and column scaling second transforms our linear system Equation (6) from the original space to the so-called scaled space, denoted by a tilde. The transformation has the form

(16)

which we abbreviate as

(17)

The notation for W stems from its Cholesky decomposition W = W ^{T/ 2} W ^1/2 (Reference WunschWunsch, 1996, p.159). Similarly, S has the Cholesky decomposition S = S ^{T/ 2} S ^1/2 and contains the column norms of the already row-scaled matrix W ^−T /² M on its diagonal.

The SVD is applied in the scaled space. Back transformation of the solution in the scaled space to the desired solution v in the original space is carried out by . It can be shown that for a full-rank underdetermined (over-determined) system, row (column) scaling is irrelevant as the respective scaling matrix is no longer present in the solution (Reference WunschWunsch, 1996, p.161, 164). Despite this fact, we always apply both scalings to cover all general cases. In addition to scaling, the use of W and S allows a degree of control of the relative norms of solution and residual.

3.5. Separation of mean and variation

Depending on the problem we are dealing with, information about the variations of the velocity around an average is more interesting than the average velocity, as the velocity variations tell us more about the processes occurring at the ice-sheet surface and their interaction with ice dynamics. Unfortunately, the minimum-norm property of the SVD will result in a solution that is smallest, in the sense of being closest to zero. This means that we might derive the wrong velocity field structure. It is therefore feasible to consider only the variations of the flow field on a homogeneous background. Hence, we separate the mean flow from its spatial variations by

(18)

where is the mean flow field and is the spatial variation. Separate mean values and , each averaged over the entire domain, are used for horizontal and vertical velocities, respectively. and , where i _n is a vector of length n with all ones. Our linear system Equation (3) can then be reformulated as

(19)

In the case that the mean velocities used for this separation are incorrect, the SVD solution of the inverse problem will try to correct this error (e.g. by providing a velocity variation on average very different from zero). For the rest of the paper we drop the tilde. We assume that separation of mean and variation and subsequent scaling has been applied prior to SVD. The results are then discussed in terms of the variational component of the velocity field , as well as the complete velocity field v.

4.1. Scenarios

Synthetic scenarios of flow are created with the forward model, using physical parameters chosen to mimic real conditions. The horizontal flow field u ^ref is prescribed. A Gaussian variation in surface accumulation is superimposed, i.e.

(20)

Where is the background accumulation, a typical value for the Antarctic plateau. The maximum accumulation occurs at x_μ = 0.5(x _min − x _max), the centre of the x domain, with determines the width of the distribution (Fig. 2a). Following Reference Richardson and HolmlundRichardson and Holmlund (1999), density is parameterized as

(21)

The variables ρ ₀ = 400 kg m ^{− 3} and ρ _i = 900 kg m ^{− 3} represent the density at the surface and the density of solid ice, respectively, and c_ρ = 0.05 m⁻¹. Such a density distribution is commonly observed in Antarctica.

Fig. 2. (a) Accumulation forcing and (b–e) resulting age–depth distributions using different horizontal velocities. Scenarios include: (b) no flow, NF; (c) slow flow, SF; (d) moderate flow, MF; and (e) moderate divergent flow, MDF, for the upper 50 m of the firn column (Table 1). Greyscale represents age values at grid nodes, with the spatial resolution of the greyscale corresponding to the resolution used fordiscretizingthe inverse problem. Contours are lines of equal age. Horizontal flow is from left to right. Crosses in (a) indicate position of nodes on A-grid and scale on the right is vertical velocity at the surface.

For the numerical forward model and the inverse problem, the continuous functions defined in Equations (20) and (21) are discretized onto the respective grids. The triplex-staggered grid used in the inverse problem of the linear system Equation (3) was explained in section 2.4, with more specifications given below. The forward model is implemented on a grid spanning 5 km in the horizontal and 100 m in the vertical direction, containing 51 × 101 nodes (Table 1). This volume is sufficient to cover the firn region of cold polar or high-altitude sites. It also comprises those length scales which show prominent variations in internal-layer architecture over short distances, as imaged by radar at various places in Antarctica (Reference Rotschky, Eisen, Wilhelms, Nixdorf and OerterRotschky and others, 2004; Reference Arcone, Spikes and HamiltonArcone and others, 2005; Reference Anschütz, Eisen, Rack and ScheinertAnschütz and others, 2006).

The effect of four different flow regimes of firn with prescribed horizontal velocity field (Table 1) on the age–depth distribution is displayed in Figure 2. In the simplest case, no horizontal advection takes place (scenario ‘no flow’, NF). This could be considered the case on a broad ice dome or along an ice divide. The other cases consider constant slow flow (SF) and constant moderate flow (MF) , which are also typical for polar ice sheets (Reference Bamber, Hardy and JoughinBamber and others, 2000; Reference Xiaolan and JezekXiaolan and Jezek, 2004) or high-altitude alpine glaciers (e.g. Reference Lüthi and FunkLüthi and Funk, 2001; Reference Schwerzmann, Funk, Blatter, Lüthi, Schwikowski and PalmerSchwerzmann and others, 2006). For these three scenarios the prescribed velocity variation .

For the moderate velocity of , a fourth scenario considers divergent flow (MDF) of the form (c_u is such that u(x) increases by 20% from to over the x domain); therefore . A scenario with non-constant horizontal velocities is the most likely case to be encountered in reality, so it will be the special focus of the analysis in section 5. Typical velocities for fast ice-stream flow are not taken into account in the main part of this feasibility study, but a set-up with a higher flow velocity of 50 m a⁻¹ is treated in the application of the inverse approach to glaciological problems in section 6. The scenarios clearly show how the varying horizontal advection affects the resulting age–depth distribution (Fig. 2). For scenario SF, the effect of the accumulation variation tapers off before an affected ice particle leaves the model domain. For both MF scenarios, advection is larger so the accumulation effect is still present at the outflow of the model boundary.

4.2. Additional constraints

A standard approach to determine the parameters of a physical model, assumed to be a compatible description of a system, is to minimize an objective function that gauges the misfit between measurements and model results. Model physics are usually enforced as constraints on the minimization in the form of exact equations, so-called hard constraints (e.g. Reference WunschWunsch, 1996). For ice-flow modelling, this was presented by Reference MacAyealMacAyeal (1993) in the case of estimating the basal friction of an ice stream and later applied to real data (Reference MacAyeal, Bindschadler and ScambosMacAyeal and others, 1995; Reference Vieli and PayneVieli and Payne, 2003; Reference Joughin, MacAyeal and TulaczykJoughin and others, 2004; Reference Larour, Rignot, Joughin and AubryLarour and others, 2005). Reference TrufferTruffer (2004) estimated the basal velocity of valley glaciers.

In addition to the basic physical description of a system, certain aspects of a solution such as structure, norm or boundary values are also sometimes known a priori. This information is valuable and helps to restrict the lack of uniqueness in solutions of inverse problems. It can be included in the objective function either as a hard constraint by Lagrange multipliers, or as a soft constraint by trade-off between the norm of the solution and the norm of the data mismatch. The trade-off can be implemented in several ways (e.g. by weighting, tapered least squares or damped least squares (Reference MenkeMenke, 1989, p. 52)). Although the SVD does not explicitly employ an objective function, constraints can likewise be imposed. An example is provided by Reference Waddington, Neumann, Koutnik, Marshall and MorseWaddington and others (2007), who also use SVD to invert a linear system of equations representing a thermomechanical ice-flow model.

Each of the different sets of constraints applied in the following exercises with a synthetic scenario can, in reality, also be determined from measured data. For the problem addressed here – the flow and deformation of firn – a first guess of the flow field at the surface is usually made. Horizontal surface velocities can be measured directly (e.g. ground-based global positioning system surveys of stakes) or indirectly (e.g. observations with satellite-based interferometric synthetic aperture radar). Here, the reference velocity field v ^ref represents possible measurements and thus provides a priori information about various velocity characteristics. (For real field applications, these v ^ref would be subject to measurement errors. For the synthetic scenario, however, they are the true values.)

It is therefore possible to prescribe the horizontal velocity at one or more positions at the surface (z = 0). From this point onwards, discrete index notation is used. The surface corresponds to index k = 1 = k ₀, so that

(22)

can be prescribed on one or more horizontal nodes i at the surface. In addition to the velocity, other properties such as the derivative of horizontal velocity (e.g. uniform, divergent or convergent flow) can also be prescribed. With denoting the horizontal difference of the horizontal reference velocity at the node (i, k) between neighbouring nodes, we can constrain

(23)

Distribution of horizontal velocities with depth is deducible from measurements of borehole deformation, enabling us to use k ≠ k ₀ in Equation (22) for values at depth at the borehole location i = i _b. We can also infer properties of surface-parallel shearing, i.e.

(24)

(25)

where is the vertical difference of horizontal reference velocity at (i _b, k). The case , i.e. constant horizontal velocity along the vertical, is commonly referred to as plug flow. This case is used in sections 5 and 6 below.

Not only can horizontal deformations be deduced from borehole deformation, it is also possible to directly determine the vertical velocities by different methods. One way is to observe the movement of markings in a borehole wall (Reference Hawley, Waddington, Lamorey and TaylorHawley and others, 2004; Reference Schwerzmann, Funk, Blatter, Lüthi, Schwikowski and PalmerSchwerzmann and others, 2006). This provides similar information for the vertical velocities, i.e.

(26)

(27)

To infer information about the properties of the problem posed here, such as stability of the solution and general solution structure, different combinations of the equations constraining the linear system Equation (3) are used to increase the degree of determinacy. The constraints are enforced by expanding the number of rows of the model matrix M and the data vector d in Equation (6). Each combination of constraints is referred to as an inverse problem, which is then applied to a simulation scenario (Table 2). The simplest case (denoted Plain) does not employ further constraints and simply considers equations for age advection and conservation of mass. Other constraints are set up by prescribing conditions for u or w: the horizontal or vertical velocity at the surface as boundary condition (denoted Bu or Bw, respectively), plug flow (Pf) and horizontal divergence (Du). Moreover, combinations of these constraints are also used in the inverse problems (e.g. BwPf, BuPf, BwDu, BuDu and BwPf).

Table 2. Prescribed constraints and system properties

The inverse problem Plain shows that the principal property of the kinematic approach is underdeterminacy, i.e. there are fewer known equations M than unknowns N (M = 162 < N = 180). All other inverse problems with constraining equations are less underdetermined, with the majority being overdetermined systems (Table 2). Only the rather complex MDF scenario (moderate flow with divergence) is solved with several constraints and is used in section 5 below to discuss the solution properties in detail.

The SVD inversion is implemented with the linear algebra package (LAPACK) routines integrated in Matlab^TM. As most of the densification of snow takes place in the upper part of the firn column, the inverse problems only address the upper 50 m. The grid used for the inverse problems spans 11 × 11 nodes, with increments of 500 m and 5 m in the horizontal and vertical, respectively. It has a five-fold lower resolution, but its nodes coincide with a subset of the grid used in the forward problem. As a result, the fields of age and density input to the inverse problems do not have to be interpolated. A linear interpolation of the u and w reference velocity fields (u ^ref and w ^ref) is carried out to project these values onto the triplex-staggered grid (Fig. 1). Evidently, the lower resolution and the interpolation will have some influence on the results. However, this effect could be considered equivalent to small measurement errors for real data. The influence of data errors on the results is considered at the end of section 5.

5.1. Singular-value spectrum

We focus on four inverse problems with different constraints to determine the solution for the velocity field of the MDF scenario: the underdetermined and simplest case Plain; the almost-determined inverse problem Bw (boundary conditions of w at surface prescribed); and the overdetermined inverse problems BwPf and BwDu (as Bw, additionally with plug flow Pf or horizontal divergence Du prescribed as constraints, respectively (Table 2)).

The first third of the ordered singular values (up to index 72 in Fig. 3) is basically identical for all inverse problems. Beyond this index, up to index 170, the spectra of the overdetermined inverse problems fall off slowly in several steps up to index 170, whereas the underdetermined inverse problems show only one or two further steps before falling steadily.

All spectra show a final discrete drop at singular values of ∼0.25–0.5. Such an abrupt and final discrete drop in a singular-value spectrum is a typical phenomenon for various problems (Reference MenkeMenke, 1989). Beyond the final discrete drop, all spectra fall continuously on the log scale. The spectra for underdetermined inverse problems decrease with increasing index faster than the overdetermined inverse problems. Whereas the rate of decrease of the spectrum for Plain does not change significantly, the other inverse problems show an increasing rate of decrease for the smallest singular values on the log scale.

In general, the spectra show greater differences for approximately the smallest 20–30% of the singular values. This has important implications for the residual norms and solution norms. Using the untruncated spectra for estimating the model parameters usually results in very small residual norms, equivalent to high parameter resolution, but larger solution norms. The corresponding velocity fields show very detailed velocity structures which, however, need not be correct.

To demonstrate the influence of the choice of the reduced rank , Figure 4 displays the resulting difference norms for the inverse problem Bw of the whole range of possible values for . The difference norm of vertical velocities , weighted with the square root of its mean, is constant at about 1 for , then falls off rapidly to steady values around 0.26 before it rapidly increases for . This distribution indicates that the vertical reference velocity structure is best approached for , although not exactly matched. This can be confirmed by checking the complete velocity structure for other in figures comparable to Figure 6 (omitted for brevity).

The distribution of for Bw, likewise weighted with the square root of its mean, is constant around 1.7–1.8 for . Two plateaux are present for . In this region, the mismatch of horizontal reference and solution velocities is at its minimum. For , increases with rank .

Adding both velocity-difference norms, each weighted with the square root of its mean velocity, a broad minimum with two plateaux for is apparent again for . This range corresponds to the tail of the singular-value spectrum (Fig. 3), where the singular values fall continuously. Similar analysis for the variation of difference norms with reduced rank for the other inverse problems yields equivalent findings: the velocity-difference norms always show a minimum for a range of singular values before showing the tendency to decrease more rapidly with larger indices. Within this minimum region, the choice of leads only to small differences of the final velocity solution. The inverse problems which constrain the horizontal velocity at the surface, i.e. Bu, BuPf and BuDu, basically display the same features.

One choice for is the index of the last step-like drop-off as the lower bound of the singular-value spectrum, used for estimating the solution of our inverse problem in Equation (9). The continuously and rapidly falling part of the singular spectra is thus truncated, a common practice when using SVD for solving inverse problems (e.g. Reference WunschWunsch, 1996).

This leads to poorer resolution, but smaller solution norms and velocity-difference norms, and yields sufficiently realistic results for most inverse problems (Fig. 6). Although the resulting smallest singular value of the truncated spectra is about the same order of magnitude for all inverse problems, the corresponding reduced rank differs significantly (Table 2). This results from the fact that, depending on the number and type of constraints, the equations show a varying degree of linear independence. The smaller the singular values, the more linearly dependent are the equations. For BwPf, however, this choice of at the final discrete drop produces a field of almost constant horizontal velocities, implying that important information is still present in the tail of small singular values for larger . For BwPf it is actually possible to maintain the full rank and obtain realistic solutions.

To accommodate this observation, a more subjective choice for would be to choose a singular value of 0.2 as the cut-off value for all inverse problems. Which choice to make is in general difficult, especially if no further details are available from a priori information. The final discrete step-like drop-off was chosen for all inverse problems, apart from BuPf and BwPf for which full rank is maintained. This is justified as BuPf does not display a falling tail of singular values (not shown) and the drop for BwPf occurs only for a very large index and is less severe than for the comparable spectrum of BwDu (Fig. 3).

5.2. Model and data resolution

The resolution matrices T_U and T_V provide other means of judging the solution of an inverse problem. If non-diagonal elements are non-zero the related main-diagonal element must be less than unity, indicating that this parameter is not fully resolved, i.e. only averages of nearby parameters can be determined.

The solution of three inverse problems with different constraints for the MDF scenario is now discussed. At full rank, the data are fully resolved for all underdetermined inverse problems, and the model parameters are fully resolved for all overdetermined inverse problems. The latter is the case for BwPf, for which the full rank R = 180 is maintained (Fig. 5). For the truncated underdetermined solutions (Plain and Bu) discussed for the MDF scenario above, the model resolution matrix T_V indicates that the horizontal velocities are only poorly resolved (Fig. 5b). The vertical velocities are equally well resolved for both inverse problems. It will become evident that this is in accordance with a comparison of the actual velocity fields shown in Figure 6. Without checking the reference velocity field it is therefore possible to assume that, in the underdetermined cases, the vertical velocity solutions are more reliable than the horizontal velocity solutions.

For all overdetermined or truncated underdetermined cases, the data cannot be fitted exactly. This gives rise to larger residuals. The order of the diagonal elements of the data-resolution matrix T_U in Figure 5a follows from that of the structure of M of the linear system Equation (6), rearranged as a vector. This vector represents groups of equations or constraints (as indicated on the abscissa of Fig. 5a). Within each group of equations, the elements are sequentially ordered by horizontal rows of grid nodes. The diagonal elements now indicate that the data predicted by the age equations are poorly resolved for all inverse problems.

The equations for conservation of mass are better resolved, though not fully. In particular, they show a decreasing resolution trend with depth (larger element index) without further constraints (Plain). For BwPf, the plug-flow constraint is well resolved. This result is evident in the horizontal velocity structure discussed in section 5.3. For both Bw and BwPf, the equations prescribing the vertical velocities at the surface are very badly resolved. The oscillations in data resolution are not arbitrary. The variations visible in Figure 5a seem to systematically depend on the position of the underlying node. The variations are smaller in the horizontal than in the vertical direction. Overall, the model-parameter and data-resolution matrices allow us to judge and improve the quality of the solution by inspecting the residual and solution norms and the singular-value spectrum without requiring a reference velocity field.

5.3. Solution vs reference velocity fields

The principal results obtained in section 5.2 are clearly seen in the velocity distribution (Fig. 6). The reference velocity fields u ^ref and w ^ref, which are the correct solutions being sought, are shown in Figure 6a and a’. The underdetermined solution Plain without constraints does not reproduce the horizontal velocity, but gives an idea of what the vertical velocity field might look like. In the almost-determined case Bw, the vertical structure is reproduced correctly but the vertical velocities in the solution are smaller than the reference velocities. The horizontal velocities again do not show the expected divergence. The vertical velocities in the overdetermined case BwPf are very similar to the almost-determined case Bw, but differ slightly more from the reference values.

Because plug flow was used as a constraint for this inverse problem, the horizontal velocities u ^est are now in very good agreement with the reference field u ^ref, although with smaller values overall. The better agreement of the horizontal velocities of the solution and reference is consistent with the fact that the horizontal velocities are well resolved for this inverse problem (see diagonal elements of model resolution matrix in Fig. 5b).

5.4. Norm properties of solutions

We next discuss the different norm properties of the inverse problem with different constraints as listed in Table 2. The difference norm for horizontal velocities is very sensitive to the choice of the mean velocity . To provide a similar foundation for all inverse problems, the mean velocity is always provided as the mean of the reference velocity field for each scenario, such that only the variations in the velocity solutions are compared (Table 1). The influence of zero-mean velocities is discussed below.

For full-rank SVD, ordering the different inverse problems with increasing M (the number of equations) as in Figure 7 would generally illustrate the dependence of the residual norm on determinacy. Naturally, for full-rank underdetermined systems (M < N) the data can be fit exactly, resulting in . For reduced rank, however, the residual norm increases, but yields a smaller solution norm . For the MDF scenario, the residual norm is more than a factor of two times larger for the underdetermined problems (Plain, Bw, Bu) than for the overdetermined problems (BwPf, BwDu, BuPf, BwDu) (Fig. 7a). In each of these two groups, the residual norm is quite constant. The velocity norm spans an order of magnitude (Fig. 7b), with the opposite ratio for under- and overdetermined inverse problems than for the residual norm, as expected.

More interesting from an application point of view is the residual between reference and solution velocities (Fig. 7c and d). The difference norms of horizontal velocities drop from for underdetermined problems to values close to zero for overdetermined problems. The difference norms of vertical velocities vary around for the over- and underdetermined problems, except for the cases Plain and BuDu which are only slightly larger with .

In some experiments, a priori information on horizontal velocity fields may be unavailable. In those cases, would have to be used. Employing this case for the MDF scenario, the velocity-difference norm remains quasi-constant, but the residual norm significantly increases for those inverse problems that do not incorporate boundary values for u at the surface. Without a non-zero estimate for mean velocities, the solution produces the smallest velocity norm as a consequence of the minimum-norm property of the SVD. Reducing the rank does not provide a remedy in this case.

5.5. Error and covariance estimates

The final point to investigate, fundamental to all inverse problems, is the solution uncertainty. The quantities density ρ and age A are part of the data kernel M. Density measurements along ice cores are very accurate, usually with an uncertainty of <2%. However, our assumption of a laterally homogeneous density distribution might be wrong, even if mean distributions are considered. The uncertainty of the age–depth distribution determined from radar surveys depends on factors including: conversion of radar travel time to depth based on integrated density; estimating age from ice cores; transferring the ice-core age to the internal horizons; tracking of individual horizons; and interpolation of the age distribution onto the SVD grid.

From analysis of Antarctic field data, Reference Eisen, Nixdorf, Wilhelms and MillerEisen and others (2004) found a maximum error of approximately 2% for the age–depth distribution in firn. In alpine regions, or regions with a laterally inhomogeneous density distribution, this error might be larger.

An error estimate of the model parameters requires knowledge of the data covariance C_vv , according to Equations (12) and (14). For the linear system considered here, uncorrelated uniform variance for the data cannot be assumed, as different physical equations are taken into account. Instead of prescribing an arbitrary data covariance we perform a Monte Carlo-based estimate of covariances, using perturbed reference velocities, age and density distributions as input to a forward calculation using Equation (6).

A total of 10³ experiments, each of which uses a normally distributed random error of 10% for A, 2% for ρ and 1% for v ^ref, results in a distribution of estimated data vectors. From this, the corresponding distribution of residuals n follows. Subsequent analysis finally yields an estimate of the residual covariance R_nn . As could be expected from the numerical set-up, the different equations are not uncorrelated. Although the main diagonal dominates R_nn , secondary diagonals also exhibit significant components. The contribution of the covariance of the nullspace vectors through R _αα to the model uncertainty is neglected, as a priori information about its structure is not available.

We compare the model uncertainty for the solution obtained with inverse problems BuPf and BwPf. Following Equations (18) and (19), both inverse problems with constraints for the MDF scenario solve the velocity variation on a background velocity of and . The model uncertainties P_vv for and of the solution of BuPf at full rank increase with element, i.e. depth (Fig. 8a and b). (Note that according to Equation (5), the elements of the vector are sequentially ordered by horizontal rows of grid nodes.) For the horizontal velocity variation , maximum uncertainties occur at larger depth and are approximately equal to the maximum velocity variation (Fig. 8a). For the vertical velocity variation , uncertainties for near-surface nodes are an order of magnitude smaller than the velocity variation, and at larger depths are approximately equal to the maximum variation (Fig. 8b). The uncertainty estimates for are always at least one order of magnitude larger than the actual residual between the solution of BuPf and the reference velocity. For , residuals and uncertainties are of comparable magnitude.

Fig. 8. Elements of the solution vector and the reference for velocity variation, the solution and reference residual vector (and solution uncertainty) and the diagonal of for the inverse problems (a, a′) BuPf and (b, b′) BwPf at full rank for the MDF scenario (Table 1). The display is split into (a, b) horizontal velocities and (a′, b′) vertical velocities . In (b) the y axis on the right corresponds to the elements of P_vv , as they are two orders of magnitude larger than the velocity variation . The components of each vector correspond to sequentially ordered horizontal rows of grid nodes. For instance, the uppermost horizontal velocities of the solution domain correspond to elements 1–11 and the nodes in the row below to elements 12–22.

For the inverse problem BwPf, also solved at full rank , the uncertainty of the horizontal velocities is about two orders of magnitude larger than the maximum velocity variation (Fig. 8c). The uncertainty for the vertical velocity variation is comparable to those of BuPf (Fig. 8d).

Although BuPf and BwPf produce very similar solutions for the velocity field, the uncertainties of their horizontal velocities of the solution are very different. This can be attributed to the different constraints for the horizontal velocity. For BuPf the horizontal surface velocities are prescribed as a constraint. By constraining plug flow, the horizontal velocities at larger depth are also a constraint. For BwPf, merely plug flow is a constraint. The actual value of the horizontal velocities is therefore more influenced by the age–depth field for BwPf than for BuPf, and thus subject to larger uncertainties.

The uncertainty of the residuals n, and therefore of the model covariance C, depends significantly on the rank chosen. Generally, for close to the full rank R, the uncertainties of the solution are larger than for smaller . For example, for the uncertainties for the horizontal velocities of BwPf are comparable with those of BuPf for full rank with . Reducing the rank used for the solution leads to smaller uncertainties, but decreases the resolution of the model parameters. Again, this is the manifestation of the trade-off between resolution and model covariance. Moreover, for the covariance of the nullspace vectors R _αα contributes to the model uncertainty of Equation (14), but cannot be estimated without a priori information.

6.1. Variation of the migration velocity

Under certain conditions, an accumulation pattern is migrating at a different velocity from the underlying ice. This is the case for megadunes on the Antarctic plateau (Reference Fahnestock, Scambos, Shuman, Arthern, Winebrenner and KwokFahnestock and others, 2000; Reference Frezzotti, Gandolfi and UrbiniFrezzotti and others, 2002) and smaller dune-like features in coastal areas (Reference Anschütz, Eisen, Rack and ScheinertAnschütz and others, 2006). Although estimates of the horizontal velocity of the ice might be available, we cannot use it to deduce the migration velocity of the accumulation pattern. The internal-layer structure provides the key to the answer, however, as it is influenced by the relative velocity between the accumulation pattern and the ice and not by the absolute velocity of the ice itself. Using the surface ice-velocity constraints under such conditions would not result in a realistic pattern of vertical velocity and accumulation. It would be more reasonable to prescribe additional flow conditions and determine the migration velocity of the accumulation pattern relative to the ice surface by solving the resulting inverse problem.

Assume that from a field survey, ground-penetrating radar (GPR) data and dated firn or shallow ice cores are available. The age–depth distribution results from merging the GPR profile with the age and density profiles of the core. For brevity, let us assume that the true distribution of the migration velocity and other physical properties corresponds to the MDF scenario as treated before. We now take the shallowest internal layer as a proxy for surface accumulation and use it as the first constraint Bw. Although this internal layer is subject to advection relative to the accumulation pattern, as are the deeper layers, the advected distance will in general be small enough to provide a first guess of the surface accumulation. As we only have shallow GPR data covering the firn column, we can assume plug flow in the firn and use this as the second constraint Pf. We therefore have the inverse problem with constraints BwPf, different properties of which have been determined and discussed for the MDF scenario in section 5.5. The horizontal and vertical velocity fields of the solution to the inverse problem are those shown in Figure 6.

For the glaciological problem assumed here, the horizontal velocity field now corresponds to the relative horizontal migration velocity of the accumulation field with respect to the ice. If the ice velocity is also available, the absolute migration velocity of the accumulation pattern can be calculated. This example shows how useful a kinematic inverse approach can be in providing an estimation of the horizontal advection velocity field, even if no further velocity information is available.

6.2. Estimation and stationarity of an accumulation pattern over time

The reader may wonder why it is actually necessary to use a mathematically complex inversion scheme under the simplifying assumption of plug flow in firn. If the flow is indeed plug flow, then all information on the horizontal field could be deduced from measurements at the surface. However, determination of accumulation from the age distribution produces significantly different results for conventional techniques and for inverse-problem solutions. With the conventional technique, accumulation is estimated as the quotient of cumulative mass difference and age difference between two isochronous layers. The effect of advection on layer architecture for an inhomogeneous accumulation pattern can lead to non-intuitive results, as demonstrated for a number of cases by Reference Arcone, Spikes and HamiltonArcone and others (2005).

A spatially varying accumulation pattern and strong advection cause convolution of surface signals to appear in the internal-layer structure. From the conventional approach, even if a correction for advection is included, it is impossible to tell whether the accumulation pattern and value was constant in the past. Misinterpretations of internal-layer data are therefore possible.

To demonstrate the capability of the kinematic inverse approach in providing the answer for this case, the problem of a spatially oscillating accumulation pattern and a constant advection velocity (Reference Arcone, Spikes and HamiltonArcone and others, 2005, fig. 10c) is now discussed. This problem is comparable to the MF scenario (Table 1) with a higher horizontal velocity. For the analysis, the age–depth field (Fig. 9a) is produced by the forward model for a model domain of 30 km in the x direction and 100 m in the z direction. Cyclic boundary conditions are used at the inflow of the model domain, mimicking an infinite extension of the accumulation pattern at the surface. Accumulation and flow parameters are comparable to the scenario of Reference Arcone, Spikes and HamiltonArcone and others (2005, fig. 10c): a constant horizontal velocity u = 50 m a⁻¹, a stationary cosine-like accumulation pattern with a wavelength of 10 km, a mean accumulation of and a spatial amplitude of (Fig. 9b). The density–depth function is the same as before. We consider only the first 25 km as the model domain of the inverse problem.

Fig. 9. SVD solution vs conventional accumulation estimates and prescribed values. (a) Age–depth distribution according to the scenario presented by Reference Arcone, Spikes and HamiltonArcone and others (2005, fig. 10c) with ice flow u = 50 m a⁻¹ from left to right (note the almost horizontal isochrones for an age of around 200 years); (b) prescribed surface accumulation (black line and grey crosses) producing the age–depth distribution of (a); (c) accumulation solution for the inverse problem BuPf calculated from the vertical velocity solution and the prescribed density–depth distribution; and (d) conventional accumulation estimates with correction for horizontal advection according to layer age. Grey crosses in (b–d) indicate reference values for accumulation at numerical nodes at the surface.

For the conventional accumulation estimate, advection can simply be taken into account by shifting the accumulation distribution determined from neighbouring internal layers upstream by the distance covered with the mean horizontal flow velocity, since the layers were deposited at the surface (Fig. 9d). The result shows that, apart from the accumulation pattern derived from the layer closest to the surface, the accumulation values calculated from deeper layers vary considerably from the actual accumulation pattern at the surface.

For accumulation minima at the surface, the accumulation derived from the deepest layers determined with the conventional approach is up to 70 kg m ^{− 2} a⁻¹ higher than the actual accumulation at the surface, equivalent to ∼70% of the reference value. For accumulation maxima, it varies by about ±30 kg m ^{− 2} a⁻¹ (8.5% of the reference value). In addition, the conventional accumulation pattern cannot be reconstructed over the complete x domain, because the layer architecture essentially necessary for a complete reconstruction has been partly advected outside the domain of the known age–depth distribution (the deepest continuous layer has an age of about 340 years, corresponding to an advection of 17 km). For a detailed analytical discussion on the related topic of causal relations between changes in accumulation, layer architecture and particle trajectories, see Reference Parrenin and HindmarshParrenin and Hindmarsh (2007).

For solving the kinematic inverse problem, we assume that the horizontal surface velocities are known from measurements and that plug flow prevails. We can then use the constraints BuPf. To provide enough numerical nodes per wavelength of the accumulation pattern, it is necessary to increase the resolution of the grid to 25 × 25 nodes. This yields a total of M = 1609 equations and N = 1104 unknowns. As for the comparable case mentioned earlier, the full rank R = 1104 is applicable to the inverse problem BuPf.

The accumulation can be determined from the vertical velocities of the solution. It provides a congruent distribution for the accumulation derived from the vertical velocities at all depths. For accumulation maxima, the solution is about 5 kg m ^{− 2} a⁻¹ smaller than the actual accumulation pattern, equivalent to −1.4% of the reference values. For minima, it is about 5 kg m ^{− 2} a⁻¹ larger, equivalent to +5.2%. The congruent shape of the accumulation pattern derived by the inverse approach implies that the assumption of steady state is correct. This in turn tells us that the accumulation was constant over time. Again, this result cannot be achieved from the conventional accumulation estimates alone.