Introduction

The classical and simplest studies which are done on a temperate mountain glacier are to measure mass balance and surface velocity at a set of stakes. Since the rheological law for permanent (tertiary) steady creep of glacier ice is today more or less well known Reference Duval(Duval, 1981), we can tackle next the inverse problem, that is the calculation from these surface data of the velocities at depth. Then, in particular, if the bedrock depth has been determined by seismic exploration or boring, sliding velocities may be determined.

This problem is a very difficult three-dimensional one. As a first step only the two-dimensional problem, when the valley glacier is cylindrical and all the velocities are parallel, will be considered. Assuming an isotropic viscous rheological law (i.e. at a given point deviatoric stresses and strain-rates are proportional), the generalized Navier-Stokes equations hold. In this two-dimensional case, only one of them has to be considered. Mass balances are assumed to be negligible and the glacier surface to be a plane xOy with slope tan α. With Ox in the direction of the velocity and Oz pointing downwards. This Navier-Stokes equation is (Reference LliboutryLliboutry, 1971, equation 65)

(1)

Assuming Glen’s rheological law, the viscosity η is given by

(2)

The Navier-Stokes equation then becomes

(3)

In our calculations the value n = 3 has been adopted.

The direct problem has been solved numerically by Reference NyeNye (1965) assuming no sliding, or a uniform sliding velocity. It has been solved by Reference ReynaudReynaud (1973) assuming a solid friction law with a uniform friction coefficient and assuming the presence of water with a single well-defined piezometric surface at the ice bedrock interface. Although Reynaud’s assumption seems more realistic, it is only a rough approximation. It involves two unknown parameters and only one of them can be deduced from global equilibrium. The computation needs a large computer, and to solve the inverse problem by trial and error would waste a great deal of computer time.

Let us then proceed directly and try to solve Equation (3) numerically from the surface downwards with the initial boundary conditions:

(4)

If we can determine the velocity at each successive step in z, the computation will be simple enough to be done on a mere desk computer. It turns out that such an extrapolation downwards is unstable, contrary to what might be thought from the existence of Reference SomiglianaSomigliana’s (1925) solution in the Newtonian viscous case. Nevertheless a more careful procedure allows significant velocities to be computed down to some depth.

Two-dimensional explicit finite difference method

As a first trial, partial derivatives at the point (y = mh, z = pk), where the velocity is u_m ^p , were approximated to the second order by the following expressions:

(5)

The value u_m ^p+1 does not appear in the approximation of ∂u/∂z in order to keep the finite-difference equation linear for the unknown values u_m-1 ^p+1 , u_m ^p+1 , u_m+1 ^p+1 , at depth z = (p+1)k. The same assymetric formulae were used for the mesh-points on the left and right boundaries. For the pth line of the grid, a linear system was solved and gave the values u_m ^p+1 at the (p +1) th line. This method was programmed on an IBM 5100 desk computer and tested for the cases of an infinitely deep channel and of a channel with a semicircular cross-section. For any value of meshes h and k, a numerical instability appeared at shallow depth.

Such an instability could have been predicted using the heuristic principle of fixed coefficients, which allows von Neumann’s stability condition Reference Godunov, Rabiyenki and Moscou(Godunov and Rabiyenki, 1977, chapter 8) to be applied. Putting constants in place of the coefficients of second-order partial derivatives, and k/h = r, the homogeneous part of our finite-difference equation (which controls the behaviour at infinite distance) becomes:

(6)

It is obvious that this recurrence scheme is not stable against a perturbation if u_m ^p+1 is equal to λu_m ^p , where λ is a constant independent from p, and if │λ│ > 1. With a perturbation proportional to exp(iβm) it is easy to compute the λ(β) for which u_m ^p+1 = λu_m ^p (i.e. the eigenvalues of the operator which transforms u_m ^p into u_m ^p+1 ). Von Neumann’s necessary (but not sufficient) condition for stability is then that a full spectrum of eigenvalues exists with moduli all less than 1, for any β.

Putting u_m ^p = λ^p exp(i βm), Equation (6) becomes:

(7)

For any β there are two eigenvalues λ_1, λ₂ for which │λ₁│.│λ₂│= 1. Thus the equation with fixed coefficients does not satisfy von Neumann’s condition.

Some considerations concerning somigliana’s solution

Assuming Newtonian viscosity (n = 1), Equation (3) reduces to the Poisson equation

(8)

The corresponding explicit finite-difference equation is then precisely Equation (6) with c = 0 and a = b. Thus for a Newtonian viscosity the numerical instability of this finite-difference equation is mathematically proved. (In the general case we had only an heuristic argument in favour.) This statement will surprise many glaciologists since it is well known that Reference SomiglianaSomigliana (1925) has given the following analytical solution for the inverse problem.

Letu = u _s(y) be an analytical expression of the surface velocity. The function of a complex variable u_s (y+iz) may then be defined. Let us consider its real and imaginary parts

(9)

and consider the solution

(10)

For z = 0,u_s (y+iz) is real and then from Equation (9) u_s (y) = U(y,0). Also from Equation (9)

(11)

Thus ∂U/∂z vanishes for z = 0, and the second boundary condition at the surface η du/dz = 0 is also satisfied. The inverse problem has then been solved, theoretically at least.

Actually u_s(y) is known only at a set of stakes. Assuming that the distance between these is h, the analytical expression which is adopted for u_s (y) may be wrong by terms in sin (nπy/h), where n is an integer. Then our analytical solution may be wrong by terms in

(12)

where n is any integer. These terms would grow exponentially if z > h/(nπ). The lack of information between stakes introduces high-frequency noise, which may be removed by smoothing. Whether this optimistic conclusion remains true for non-Newtonian viscosity is an open question.

The fact that we have only a limited number of surface velocities, because the glacier width is finite is much more troublesome. The adopted analytical u_s (y) may be wrong by a polynomial of any degree n ⩾ = N — 1, The corresponding error in the solution has terms in zⁿ (if n is an even integer) or z^{n -1} (if n is odd). The same kind of error may come from any small inaccuracy in the velocity measurements (and, when using a computer, from truncation errors). Let the glacier extend between y = -L and y = +L, and ±∊ be the inaccuracy. Then we cannot assert that an error:

(13)

has not been introduced. This will lead to an error in the solution which, for y = 0, if n is even, amounts to:

(14)

This error cannot be suppressed by any smoothing procedure and means we have little chance of obtaining valuable velocities for z > L. Nevertheless, coming back to the non-Newtonian case, it would be sufficient if the computation remained more or less stable down to z = L, since the maximum depth of valley glaciers is seldom larger than half their width.

“One-dimensional” finite-difference method

In the non-Newtonian case several finite-difference schemes have been tried. The one which remains stable at largest depths is as follows.

The N given surface velocities u_s (y) are represented by a polynomial of degree N-1. In order to smooth as much as possible without loosing the input of information, velocities at any depth will also be represented by polynomials of degree N— 1:

(15)

Horizontal derivatives are then calculated directly as polynomials in y, while vertical ones are approximated by finite-difference expressions.

At the start the N coefficients a_m (0) are known and a_m'(0) = 0. Derivatives ∂u/∂y, ∂u/∂z, ∂²u/∂y² and ∂²u/∂y ∂z can then be calculated for z = 0 at the N selected values of y (say y_i ). Equation (3) then gives the second-order vertical derivative at the N points (y, 0 ).

Interpolating by a polynomial in y of degree N —1 gives the a_m" (0), since:

(16)

Let us increase z by k. The a_m (k) and a_m’(k) may be obtained by extrapolation. Equation (3) gives N values (∂²u/∂z²) y_t.k whence the a_m"(k), and so on.

Extrapolation formulae of increasing accuracy are used at first; after the third line the following equations, which are valid to order 6 and 5 respectively, are kept:

(17)

For an infinitely deep channel of width 2L, the exact analytical solution is

In this case our procedure gives correct numerical results down to a depth more or less equal to the glacier width 2L.

Fig. 1. Semicircular cylindrical channel, no sliding. Velocity versus transverse coordinate y at different depths drawn from the analytical solution (dashes) and from numerical computation (solid line).

For a channel of semicircular cross-section with radius R, an exact analytical solution is:

In this case our procedure gives the results plotted on Figure 1. Solid lines represent the numerically computed velocity profiles at different depths u_z (y), and dashed lines the corresponding true profiles. R = 450 m was adopted. For z = 330 m the numerical solution becomes meaningless. For z = 360 m numerical instability becomes so important that the computed curve could not be contained in the figure.

Lastly the computation was made using stake velocities measured on Glacier d’Argentière, where the width 2L = 900 m. Numerical instabilities appeared long before the maximal depth (350 m) had been reached. The calculated velocities seem to be correct only down to a depth of 150-200 m (Fig. 2).

Fig. 2. Transverse section of Glacier d’Argentière at about 2650 m a.s.l. Surface velocities have been measured at the plotted stakes. For numerical computation of velocities a horizontal straight line has been substituted for the actual surface. Lines within the glacier are lines of equal velocity as computed before instability happens. Velocities in metres per year.

Conclusion

Our heuristic arguments and our numerical computations allow us to assert that the inverse problem can only be solved, in the two-dimensional case, for shallow glaciers. In order to determine sliding velocities in valley glaciers, inclinometric surveys of bore holes seem unavoidable, at least until a good friction law and a good model for water pressures at the interface have been established.