1 Introduction

The effect of longitudinal strain on the state of stress at the base of a glacier or ice sheet has been studied previously by Reference LliboutryLliboutry (1958), Reference RobinRobin (1967), Budd (Reference Budd1968, Reference Budd1970) and has found its theoretical basis in the articles of Reference CollinsCollins (1968) and Reference NyeNye (1969). In all these articles attention is focused on the significance of the longitudinal stresses set up by the flow of the ice over protuberances on the bed. The idea is to improve the formula for basal shear stress which to first order is independent of the material behaviour of the ice and involves only the local glacier depth and the local inclination of its surface. The improved shear-stress formula should account for the creep deformation of the ice and hence should depend on the creep law under consideration. This is indeed the case.

According to the above result, then, local fluctuations in the variation of the upper and lower glacier surface manifest themselves in the basal shear stress (and other quantities). However, since local fluctuations cannot govern the overall mechanical behaviour of the glacier, only quantities which are averages over certain respective lengths are of physical significance. Dependent upon the size of this length, certain details in the local fluctuations are filtered out, and it is conceivable that for sufficiently extended averages the resulting mean stresses will agree with those obtained if the ice sheet is replaced by a slab of constant thickness. Such calculations were performed by Budd in his application to the Wilkes Ice Cap. His findings are that the basal shear stress does not fluctuate in sympathy with the surface slope over distances of roughly four times the mean glacier thickness, but it does follow the surface slope changes over distances of the order of 20 times the glacier depth. Consequently, the basal shear stress fluctuates much less than the surface slope.

In a later paper (Reference BuddBudd, 1971) on the distribution of stress and velocity, an analysis was performed following the theoretical investigations of Collins and Nye in order to find the range of wavelengths which would sufficiently influence the basal shear stress formula and hence would corroborate the above results. Budd’s analysis rested on his earlier paper (Reference BuddBudd, 1970) which was critically reviewed by Reference Hutter, Hutter, Legerer and SpringHutter and others (in press) and was found inaccurate. A quantification of the longitudinal strain effect therefore still remains to be given.

In Collins’ and Nye’s analyses the variation of the basal shear stress is expressed in terms of integrals of normal-stress differences and of shear stresses. These terms are unknown as long as the stresses remain undetermined throughout the glacier. Therefore, it would be advantageous if a computational approach could be given by which longitudinal strain effects evolved by means of a method of successive approximation. I introduced such a method (Reference HutterHutter, 1980) in an article on time-dependent surface elevation of an ice slope. Starting with a long-wavelength assumption, a stretching of coordinates and a subsequent perturbation approach were introduced to solve the respective equations in a step-by-step manner. Here, I use the same approach again and shall demonstrate that the well-known basal shear-stress formula may be derived as a first-order approximation of this perturbation scheme, but more importantly, I derive explicit higher-order formulae with the aid of which the effect of material properties on the stress formulae can be determined. The method has the advantage of being explicit and that variations of ice properties with depth can straightforwardly be taken into account; it involves simple quadratures, which can be performed explicitly and analytically if the material properties of the ice do not vary with position. It demonstrates the special features of Glen’s flow law and demonstrates that the creep behaviour of ice at low stretchings should be governed by another flow law than Glen’s.

Part of the reason for writing this article was to complement the calculations of Reference Hutter, Hutter, Legerer and SpringHutter and others (in press) and to give further justification that the use of stretched coordinates is appropriate and legitimate. In Reference HutterHutter (1980), kinematic wave theory evolved as a first-order approximation of this singular perturbation method. Here, we demonstrate that the well-known basal shear-stress formula can equally be obtained as a first-order approximation. Furthermore, both kinematic wave theory and the shear-stress formula can be improved using our method in an explicit manner. This obviously gives clear indications about the validity of the classical formula. Such facts should serve as a reason in favour of, and not against, the method of stretched coordinates.

Fig. 1. Geometry and definitions for the ice sheet.

2 Governing equations

Consider the slow flow of a viscous medium down a slope with slowly varying bottom surface and slowly varying thickness. Assume the medium to obey a non-linear constitutive relationship such as Glen’s flow law or any generalization of it. Let, further, (x, y) be a Cartesian coordinate system; x is down and y normal to the plane x = 0 (see Fig. 1), and assume the deformation to take place in the (x, y) plane only. The governing field equations describing plane strain and written down in dimensionless coordinates are

(1)

where

(2)

In these equations x and y are dimensionless Cartesian coordinates, (σ _x , σ _y , τ ) dimensionless normal and shear stresses, respectively, and (u, v) are dimensionless velocities. The function characterizes the type of Maxwellian fluid adopted here. For a ≠ 0 the emerging equations describe the flow of a fluid obeying Glen’s flow law, a ≠ 0 gives a slight generalization of that.

Incidentally, in the above the non-dimensionalization is the same as in Reference HutterHutter (1980); consequently, lengths, stresses, and time are non-dimensionalized with a representative thickness of the ice sheet D, with σ ₀ = ρgD, where ρ is the mass density of ice and g the gravity constant, and with the time where A, n are phenomenological coefficients of the generalized Glen flow law

Here, D_ij is stretching, t_ij ' stress deviator, t_II'² its second invariant, and a ⁻¹ a linear viscosity (a = a/(ρgD) ⁿ⁻¹). Such a non-dimensionalization guarantees that u, v, and σ _x , σ _y , τ are O(I) as long as γ is finite. It should also be noted that D is the only quantity in the above scaling provided by field observations. All other quantities follow from material properties of ice.

In ice sheets, the material properties strongly vary with temperature and hence with y. In this case we may assume n and A to be functions of the depth variable y; y₀ is a reference depth, usually the surface or the base. More generally the variations in surface and bottom geometry also give rise to an x-dependence of the temperature field. This variation is neglected here.

Equations (1) and (2) are to be solved subject to the boundary conditions, at the base y = γ _B(x):

(3)

at the top surface, y = γ(x, t):

(4)

ϕ and m are constants describing the sliding mechanism at the base, is a dimensionless accumulation rate, and α and β are the angle of bottom and surface inclination relative to the mean inclination γ. They are defind as

(5)

Henceforth, approximate solutions to the above boundary-value problem will be sought which are based on a perturbation expansion procedure. The full determination of such solutions includes the fields and the geometry of the top surface as functions of time. The perturbation technique presented below, however, allows the separation of the two. In particular, from the assumption of the bottom and top topography, the stress and velocity fields are determinable, provided, of course, that the temperature distribution is also known. This makes the solution technique especially valuable, if field records of geometry and temperature are known. It should be stated, however, that the equation for the top-surface geometry could also be obtained; this was demonstrated by Reference HutterHutter (1980) and will not be repeated here. In the following it will therefore be tacitly assumed that the surface geometry is known from field measurements. Subsequent calculations may therefore be viewed as a recipe to obtain values for stresses and velocities at given geometry.

3 Stretching transformation

It is an experimentally established fact that surface undulations in glaciers and ice sheets are slowly varying. The same is also true for the bottom topography. It is therefore only natural to incorporate this slowness assumption into the governing equations by means of a stretching of coordinates. Thus let

(6)

μ is called the “stretching parameter” and may be defined as the ratio of a mean glacier thickness to a typical wavelength of surface undulations. Another physically reasonable definition of this parameter would be the ratio of the mean thickness to the distance over which the glacier thickness changes appreciably.

Incorporating Equations (6) into the governing equations, we obtain the following: field equations:

(7)

boundary conditions:

(8)

The stretching parameter has now entered the governing equations and since I choose μ such that γ(ξ, t) has derivatives O(I) this together with the slowly varying nature of γ with x implies that μ≪1. Values for μ suggested by field data are μ=10⁻¹ to 10⁻². This fact suggests searching for solutions in powers of μ Only after the introduction of such a perturbation expansion is the proposed solution technique approximate.

The stretching transformations (6) are not new; they have indeed been widely used in other branches of fluid mechanics. Reference FriedrichsFriedrichs (1948) and Reference KellerKeller (1948) have introduced them to derive the shallow water equations in the theory of water waves, see also Reference StokerStoker (1957). A summary of research on water surface waves based on these transformations is given in the Handbuch der Physik by Reference Wehausen, Laitone and FlüggeWehausen and Laitone (1960). The same transformations have also been used by Reference BenneyBenney (1966) in a theory of thin liquid films. In all these cases the stretching transformations (6) are used as a means of deriving a systematic way of approximating field equations which were previously derived by more ad hoc methods. In these the underlying assumptions are often unclear or hidden. To the non-specialist the emerging equations then look more rigorous than when derived with a systematic approach. It is because of his unawareness of this that he will ask for a justification of the transformations (6). In fact, there is nothing that needs to be justified with the stretchings (6) as long as the emerging equations are treated with full rigour. When approximate solution techniques are used, however, one may justifiably ask about the limitations of the results, but this is no different from other approximate-solution techniques. As regards singular perturbations, results are believed to be correct as long as they “look reasonable”, see Reference Van DykeVan Dyke (1975).

4 Longitudinal strain effects

Consider the boundary-value problem given by Equations (7) and (8) under the accumulation rate =(ξ, t). Let us construct solutions to it by the power-series expansions

(9)

Substituting Equation (9) into Equations (7) and (8), respectively, expanding the resulting equations in powers of μ and collecting terms of like powers of μ, a hierarchy of boundary-value problems is obtained. The zeroth-order equations read

(10)

and must satisfy the boundary conditions

(11)

Note that the kinematic boundary condition must not be exploited here because surface geometries are assumed known.

Straightforward integration yields

(12)

where ’ denotes the derivatives of with respect to its argument and where H = γ – γ _B.

The first-order equations are

(13)

subject to the boundary conditions

(14)

Integration then reveals that

(15)

Equations (12) and (15) furnish explicit expressions for the stresses and the velocities in a glacier with slowly varying bottom and upper surfaces. Once the latter is known, the stresses and velocities follow by mere substitution.

Higher-order terms are rather cumbersome to determine; yet second-order stresses may still be calculated with reasonable effort. The reader may show himself that

Integration yields

(16)

Notice that this expression becomes singular if Glen’s flow law is used, because in that case (0)=0. This is the reason for its extension as listed in Equation (7). Here a word of caution is in order. The difficulty with Glen’s flow law is primarily not a physical, but a mathematical one connected to the proposed perturbation approach. As a→0 the solution technique fails. It would, therefore, be appropriate to search for a mathematical procedure which remains regular as a→0. Such a solution technique would have to be sought using a multiple-variable expansion. With such a technique the possible accumulation of singularities that might arise with the above procedure when approaching the limit a→0 could be avoided. Despite this fact my solution technique is reasonable, since a is bounded away from zero with a value of roughly 10⁻³. Since μ ² (terms of this order are the first ones that become singular when the limit a→0 is taken) is smaller, in general, I do not expect that corrections O(μ ²) will become invalid.

With the aid of Equations (12), Equation (16) can be shown to have the form

(17)

where

and where T ₁ and T ₂ vanish when no ξ-dependence of the temperature field is assumed.

The evaluation of the functions T _i is rather tedious even though it is straightforward. In fact, it can be shown (see Appendix) that only quadratures are needed if is restricted to the form given by the second of Equations (2). These quadratures can be performed explicitly if the material properties of the ice do not vary with position. Calculations for this case have been performed in the Appendix.

Our aim in this paper is to find explicit formulae to evaluate the effect of longitudinal strain on the basal shear stresses. With the first of Equations (12) and the first of Equations (15) this goal has been achieved in a most satisfactory way, because, as can be seen from the Appendix, only quadratures are needed to evaluate the coefficients T ₁, …, T ₇ of the representation (17). These quadratures can be carried out algebraically if the material properties do not vary with depth, and the emerging model is a reasonable one for temperate glaciers. In cold ice as in the Greenland and Antarctic ice sheets, the stress-strain-rate relationship depends strongly on depth, since temperatures vary appreciably. In this case numerical integrations are necessary, but this is no disadvantage, because electronic computations are needed anyhow.

To find explicit formulae for the basal shear stress, let us combine the first equations of both Equations (12) and Equations (15). This yields

and since this can also be written as

(18)

At η = γ _Bfor small γ , α and β this becomes

(19)

To first order and for small inclinations, the basal shear stress is proportional to the local surface inclination (α+γ) and to the coresponding depth. This is the classical formula well known to glaciologists. Neither Equation (18) nor Equation (19) are apparently general enough to account also for the longitudinal strain effects. These are contained in Equation (17) and previous lower-order stress formulae and must be evaluated from

To second order in μ this yields

(20)

where T ₁, …, T ₇ have to be evaluated at η = γ _B.

Before I present additional information regarding numerical results, let me stress once more that the term containing μ ² as a common factor does account for longitudinal strain effects. Without this term and for small base and surface inclination angles, the basal stresses would vary in sympathy with surface slopes, which is contrary to observation. The term O(μ ²) therefore accounts for exactly this difference. That the relevant term is multiplied by a very small quantity μ ² does not, however, imply that the entire term is negligibly small; this would imply that the term in braces was of numerical O(I). This cannot be guaranteed, although the limit process μ→ o is mathematically sound. As an aside remark I might mention that a similar situation occurs in the derivation of the surface wave equation (see Reference HutterHutter, 1980); the derivation of that equation shows that the coefficient accounting for diffusionFootnote ^* is μ cot γ. Obviously for small μ this term might still be O(I) as, e.g. if γ = O(μ), which is the case for all ice sheets and most glaciers. The O(μ) ² term in Equation (20) is therefore not a priori negligible unless, of course, numerical calculations provide corroboration to the contrary. This is not the case as we shall see.

5 Numerical results

To obtain quantitative information regarding the order of magnitude of the longitudinal strain effect on basal shear, a temperate ice sheet was analysed whose bottom and top surfaces were varied according to the expressions

(21)

where

(22)

in which n = (2, 3). Evidently the bottom ice–bed-rock interface oscillates sinusoidally with an amplitude ε_B≪1 and a dimensionless wavelength λ = 1. Written in physical variables the second of Equations (21) reads (bars indicate physical variables)

H and l are the maximal thickness of the ice sheet and the wavelength of the bottom undulation, respectively, and μ=H/l≪1 is the ratio of the two. The oscillating part of the first of Equations (21) must be interpreted analogously. The expression in Equation (22) for γ ₁(x), on the other hand, corresponds to

(23)

where L is the half-length of the ice sheet. For brevity γ ₁ will be called the “equilibrium surface”. Straightforward comparison of Equations (22) and (23) yields k = l/L. The top surface varies therefore according to the equilibrium geometry of Equation (22); superimposed on it are sinusoidal undulations whose wavelength is the kth part of L.

The equilibrium profile, Equation (23), has one particular mathematical disadvantage, namely that it may lead to singular second derivatives at the head. Such behaviour is unrealistic, but it manifests itself in the expression for the second-order shear stresses. We shall in the following analysis therefore exclude the neighbourhood of both head and snout.

The above representations (21) may not be realistic in a real situation, but they will give an indication of the order of magnitude of the basal shear stress. For is close to the top geometry of an ice sheet and the sinusoidal term may be regarded as a particular term of a cosine Fourier representation of the difference .Varying ε, ε_B and k thus gives a reasonably realistic picture of longitudinal strain effects. We have performed such calculations; they are described below.

Of particular interest is the correction O(μ ²) of the shear stress. It is, therefore, most convenient to write

where , and to graphically display only the coefficient . It represents the influence of the longitudinal strain effects on the basal stresses. Calculations were performed for no-slip and for a Weertman-type sliding law by varying the coefficient Φ. The exponent in the generalized Glen flow law was given the values 2 and 3, and correspondingly a assumed the values 10⁻² and 10⁻³.Footnote ^* The mean inclination angle γ was varied between γ = 10⁻³ and 2 × 10⁻¹.

In a first set of calculations, periodic undulations were set aside, Є = 0, Є_B=0, and only the influence of the equilibrium geometry was analysed. For n=2 Figure 2 shows the second-order shear-stress corrections as a function of position and plotted for various inclination angles γ. Except for very small values of γ and close to the snoutFootnote ^* the second-order shear- stress corrections are nearly independent of x and only very weakly dependent on γ. Interestingly, second-order shear-stress corrections are virtually independent of the sliding law. This is so at least for realistic values of Ø≤ 10⁻². Generally, it may be stated that longitudinal strain effects are negligible except perhaps for γ<10⁻² and in the snout region, The corresponding results for n=3 show quite the contrary, see Figure 3! As is apparent here, the dependence of on γ is more pronounced than for n = 2. In particular for n = 3 and for small γ(≤10⁻², is independent of γ, whereas for n = 2 independence of γ seems to apply for large γ. Moreover, whereas for n = 2 results were insensitive to sliding, this can, strictly, not be said here, since for Ø=10⁻² an onset of deviations from the curves for Ø<=10⁻² can be observed (the results for Ø= 10⁻² are not graphically displayed.

Fig. 2. Second-order basal shear stress for the equilibrium geometry of a glacier. "The parameters in the generalized Glen flow law are n=2, a=10⁻². The dashed lines, for which the right-hand scale is applicable, represent the same as the solid lines. In this and all subsequent figures the symbols indicate the values of γ, the key being given in Figure 4a, and the abscissae are ξk.

Fig. 3. Second-order basal shear stress for the equilibrium geometry of a glacier. The parameters in the generalized Glen flow law are n=3, a=10⁻³. The dashed lines, for which the right-hand scale is applicable, represent the same curves as the solid lines.

The strong dependence of the basal shear stresses on the exponent in Glen’s flow law is surprising and it is to be expected that it will also be seen when sinusoidal undulations of the bottom and top surfaces are superimposed. This is true, yet when this superposition is performed, the parameters k (the ratio of undulation wavelength l to the glacier length L), € and €_B enter. Dependent on their values the behaviour due to equilibrium geometry may totally be overshadowed. Calculations were performed for k=(1, 0.75,0.5, 0.25) and (€, €_B)=(0.01, 0.05, 0.1). Undulation wavelengths thus varied between a quarter to the full half-length of the glacier and undulation amplitudes varied between one-hundredth to one-tenth of the maximum ice-sheet thickness.

In order to demonstrate that bedrock undulations of very small amplitudes contribute significantly to the second-order basal shear stresses, we have displayed in Figure 4a to d values of for the case when n=3, €_B= 0.01, €=0 and for four different values of k. When compared with Figure 3 it is seen that with decreasing wavelength the influence of the equilibrium geometry becomes smaller and smaller; for k=1 the behaviour as displayed in Figure 3 can still be recognized, although the bottom protuberance is clearly observable in the undulating behaviour of the stresses. In particular, the order of magnitude of the maxima of the basal shear stresses is similar. With decreasing wavelength these maxima become smaller and smaller; for k=0.25 the equilibrium behaviour is completely overshadowed by the bottom geometry. Nonetheless the variation of the surface geometry is not negligible in this case, because the curves in Figure 4d are far from periodic. It is, further, recognized that, for large γ, amplitude maxima drop with decreasing k. This is an indication of the fact that basal shear stresses vary in sympathy with bottom undulations very much more when wavelengths are large than when they are small.

The above results are valid for n=3. It is interesting to test whether an equally strong dependence on the wavelength of bottom protuberances would also be visible for n=2. For k=1 second-order basal shear stresses resemble pretty much the behaviour of Figure 2; with decreasing wavelength the significance of bottom undulations becomes more and more important. As a corroborating example Figure 5 (corresponding to Figure 4d) may serve, in which €_B= 0.01 and k=0.25. Except for γ=0.001 and then only close to the snout, surface geometry seems to be of little significance. On the other hand, a comparison of Figures 4d and 5 clearly indicates the importance of the value of the parameter n.

Fig. 4. Second-order basal shear stresses when a sinusoidal bottom protuberance is superimposed on the equilibrium surface geometry. The parameters in the generalized Glen flow law are n=3 and A=10⁻³. Bottom amplitude is €_B=0.01 and the wavelength parameter is for a, k=1; b, k=0.75; c, k=0.5; d, K=0.25

Fig. 5. Second-order basal shear stresses when the parameters in the generalized Glen flow law are n=2 and a=10⁻². Undulation amplitudes are €_B=0.01 and wavelength parameter is k=0.25.

It is not far from reality to state that the importance of bottom undulations also grows with growing amplitude €_B, Calculations were therefore also performed for €_B=(0.05, 0.1), and they corroborate the above statement. For k<1 this proof need not be demonstrated, because Figures 3 and 4 are sufficiently convincing. For k=1 Figures 6a and b may serve as examples. They show the ξ dependence of for €_B= 0.1, k=1 and n=2 and n=3, respectively. These figures are not only interesting for the fact that very little of the behaviour as shown in Figures 2 and 3 is left; they simultaneously demonstrate that critically depends on the exponent n in the generalized flow law. First, the tendency of the sign of as a function of x is different. For n=2, is negative in the middle part of the half glacier and positive in the head and snout region; for n=3, is positive close to the snout and negative elsewhere. Secondly, and very roughly, second-order basal shear-stress corrections are more than twice as large for n=3 as they are for n=2. These corrections are also the bigger the larger the mean inclination γ. Calculations with k <1 have also indicated that this difference diminishes with decreasing k, but the above-mentioned sign difference appears to be fairly persistent for 1≤k≤0.25. I regard this as serious, because it makes calculations regarding the effects of longitudinal strain on basal shear questionable unless one has firm knowledge about the material behaviour.

Fig. 6. Second-order basal shear stresses for the case when ϵ_B=0.1 and k=1. The parameters in the generalized Glen flow law are n=2, a=10⁻² (cf. Fig. 6a) and n=3, a=10 ⁻³ (cf Fig. 6b).

It is interesting to see whether or not similar conclusions would emerge when ϵ_B=0, but ϵ ≠0. Calculations were performed for ϵ=(0.01, 0.05, 0.1) and the same values of the wavelengths. These calculations indicate that the qualitative behaviour of the curves for n=2 and n=3 is similar. In particular, for both choices the sign of is the same in roughly the same regions (of the variable ξ) and, furthermore, second-order corrections are the bigger the larger γ. Yet, the size of these corrections is very much greater for n=3, than it is for n=2, especially when γ>0.05. Evidence for this is provided by Figures 7a to d which shows for ϵ=0.1, n=(2, 3), and k=(1, 0.25). Evidently, and quite opposite to the behaviour displayed in Figure 6, second-order basal shear stresses for n=2 and n=3 are “in phase”. Amplitude maxima grow with growing γ, but for large inclination angles γ and for n=3 they are bigger than for n=2 by an order of magnitude. This is just another indication that longitudinal strain effects should only be taken into account when firm knowledge about the material behaviour is available.

Fig. 7. Second-order basal shear stresses for the case when ϵ=0.1; a and b show the results when k=1, c and d those when k=0.25. The parameters in the generalized Glen flow law are n=2, a=10^-−2 (a and c) and n=3, a=10⁻³ (b and d). Notice the differences of order of magnitude for n=2 and n=3 respectively.

Before I conclude, mention should be made that the above analysis assumes the surface geometry to be known. Mathematically this is a disadvantage as the determination of the surface profile should be part of the problem to be solved. This could be done by substituting the results obtained for the velocities into the kinematic surface condition, the first of Equations (8), and solving the associated differential equation. However, this was not my intention here.