1.Introduction

The overall motion of a temperate glacier consists of a relative deformation or flow through the ice mass together with basal sliding over the bed (which may make a significant contribution to the observed surface velocity). The internal flow is governed by the constitutive response of the ice whereas basal sliding depends on conditions at the ice-bed interface. Specifically it is the presence of a thin water layer, caused by pressure melting and refreezing as the ice flows over undulations, which "lubricates" the bed and allows slip to occur.

The regelation mechanism was proposed first by Reference WeertmanWeertman (1957) who modelled the bed as a regular array of rectangular obstacles. The latent heat released by freezing on the low-pressure down-stream face is conducted through the obstacle to melt the ice on the high-pressure up-stream face. An energy argument and an assumed value for the shear stress at the bed then lead to a basal-sliding velocity due to pressure-melting which is expressed in terms of the obstacle and spacing lengths. An estimate of the creep-rate is made by using Glen's law for the ice response and by making further assumptions about the stresses on the obstacle. Next, an expression for the sliding velocity due to a creep process is inferred.

The two estimates of sliding velocity, one based on regelation and the other on creep, decrease and increase with the length scale of the obstacle, respectively. Weertman proposes, therefore, that since all length scales will exist in practice, sliding is controlled by the length giving the same velocity for both mechanisms. Equating the two expressions gives the sliding velocity in terms of the drag and a roughness parameter—a ratio of obstacle length to spacing. The theory is extended in a similar manner (Reference WeertmanWeertman, 1964) to account for cavitation on the down-stream faces. Later, Reference WeertmanWeertman (1971) defends this simple "semi-quantitative" theory against the criticisms of Reference NyeNye (1969,1970) and Reference KambKamb (1970). However, the theory of Weertman relies crucially on estimates of stress and a sliding velocity deduced from an inferred mean creep rate. There is also the difficulty of relating obstacle size and spacing to a more general bed profile; further, there is no demonstration that a flow field with the required features exists. A more elaborate treatment in the same spirit is proposed by Reference LliboutryLliboutry (1968) who considers a sinusoidal bed profile, accounts for cavitation, and discusses various basal friction laws, in other words, the relationships between the mean drag, the normal pressure, and the sliding velocity.

Reference NyeNye (1969,1970) and Reference KambKamb (1970) focus upon the need for an exact flow and heat conduction solution which satisfies the regelation conditions at the ice-bed interface. For this purpose Nye assumes that the ice response can be approximated as a Newtonian (incompressible) fluid of high viscosity and that the bed profile is periodic with small slope. Because of this the Reynolds number is very low, inertia terms are negligible, and slow steady viscous flow is assumed. Further, the water layer thickness is negligible compared with the length scale of the undulations, so the layer is treated as a surface distribution of heat sources (sinks) which coincide with the bed surface and which give rise to the latent heat of freezing (melting). The interface conditions are linearized to define a half-space problem for die ice flow after ignoring the upper glacier surface and introducing a small slope parameter. ε This problem is treated by Fourier analysis, first for a plane flow and then for three-dimensional flow. Second-order expressions in ε are obtained for the flow field over a sinusoidal bed. Expressions for the mean drag are deduced and the results are extended by the statistical analysis of a general bed profile. The solution assumes that no cavitation occurs. Reference KambKamb (1970) obtains the same linear solution, demonstrating that the ice motion may be neglected in the heat conduction, but uses this as a starting point for an approximate treatment when the ice satisfies the more realistic non-linear Glen law.

The present paper complements the solution of Nye and Kamb for Newtonian plane flow (small bed slope) by incorporating explicitly the depth h of the glacier and the inclination α of a (mean) bed line to the horizontal, both parameters assumed uniform over the length scale of interest. Thus, gravity (which is the driving force for the glacier motion) is included in the balance equations for slow viscous flow, and a flow solution is determined without the need to introduce an artificial shear stress at some distance from the bed (Reference NyeNye, 1969). Furthermore, the solution satisfies an upper-surface condition of normal atmospheric pressure p_a within the approximate expansion scheme. It is supposed that there is no net longitudinal stress gradient, only a periodic variation induced by the bed undulation, and that the flow is steady, that is, no rigid body acceleration of order g sin α parallel to the bed occurs; this is consistent with normal flow conditions. The mean drag on the bed is therefore

where ρ is the ice density.

Since shear stress in a thin water layer at the bed is negligible compared with τ, one boundary condition on the bed is that of zero tangential traction. This is an idealized condition which assumes that neither draining nor pinching-out of the water layer occur, and that there is no cavitation so that the ice boundary is everywhere the bed surface. The drag is the longitudinal resultant of the pressure acting over the undulations. The other bed conditions link the ice flow and heat conduction through the pressure-melting and regelation relations. For completeness a uniform geothermal heat flux acting normal to the bed is included at large distances from the bed, so that its effect on the flow can be estimated for levels which are observed typically. The effect is shown to be small. The flow solution determines explicitly the basal sliding velocity U_b, and surface velocity U_s, in terms of h, α, and the bed geometry for given ice and bed properties. An illustration for a sinusoidal bed profile shows that U_b is very sensitive to the wavelength of the undulation. The variations of U_b with seasonal changes of melt water (Reference LliboutryLliboutry, 1968) cannot be described by this model of the bed conditions.

The small parameter ε « 1 is chosen as the maximum bed slope (relative to the inclined bed line), and for mathematical completeness the bed profile is extended periodically to infinity. Dimensionless coordinates are introduced through a length scale λ defined so as to make the maximum amplitude of the dimensionless bed profile unity from a given bed line. The profile is assumed to be sufficiently differentiable for expansion procedures to a required order in ε The coefficients of leading terms are restricted to order unity by this choice of coordinates. it is assumed that λ/h is sufficiently small for the upper surface conditions to be satisfied to the required order in ε. in fact, the perturbation in the flow from a basic laminar flow is found to decay exponentially with height so that this is not a strong restriction. The approach here is to seek half-plane solutions to the perturbation flow equations in the ice and the heat conduction equations in the ice and bed which satisfy, to the required order in ε, the boundary conditions on the actual bed surface. That is, the bed surface conditions are not applied on the boundary of a half-plane as in the earlier solutions, which can therefore be correct only in the leading term of the expansions.

A complex variable method is used and solutions to order ε² obtained explicitly, though the expansions given determine valid solutions to order ε³ with a further iteration. The evaluations are elementary for a profile which can be adequately modelled by a low-order truncated Fourier series, and the solution is illustrated for a sinusoidal bed. A simple cavitation criterion is found which depends only on α and ε To an order of the first power of ε, the velocity-field perturbation agrees with the solution of Reference NyeNye (1969). However, whereas the present ε² terms are bounded, Nye has an unbounded term due to the artificial shear stress at infinity which is introduced to drive the motion. The determination of a parameter to order one, which in turn determines the basal sliding velocity, occurs in the balance of order ε² terms, so it is essential to construct a consistent expansion scheme to order ε². Furthermore, with classical boundary conditions of zero velocity there are no bounded half-plane solutions for flow which are non-trivial (Reference LangloisLanglois, 1964), so the present solution verifies that the present boundary conditions define a well-posed How problem.

2.Bed Profile and Boundary Conditions

Figure 1 displays the idealized flow problem. Coordinates (x,y) are chosen along a bed line and normal to the bed respectively, with the x-axis inclined at angle α to the horizontal. The glacier flows in the x-direction, U_s is the surface velocity at y = h, and (U_b is the basal sliding velocity defined as the x-velocity in a flow continued onto y = 0. The bed profile is given by

where f(x) is smooth, sufficiently differentiate for subsequent expansions, and supposed to extend periodically as (x) →| ± ∞·This allows a well-posed mathematical problem; a different profile distant from the region of interest will not affect the local flow field. f/(x) can be described adequately by a low-order truncated Fourier series the later solution is easily evaluated. Following Reference NyeNye (1969,1970) and Reference KambKamb (1970) it is assumed that the bed slope is everywhere small; thus,

Fig. 1. Glacier flow over a wavy bed.

Now we introduce the dimensionless coordinates (X, T) with a length scale λ:

in which the bed profile becomes

We then choose

Where

Thus

That is, the length scale λ is defined so that the contributions F(X) and F'(X) to the expansion coefficients are necessarily of order unity. The only restriction on the amplitude f_m is Equation (6) with the strong inequality represented by Equation (3) . My approach is to seek flow solutions in the half-plane y > 0 (Y > 0) which give correct values on the surface Y = Y_b, to a required order in ε. Strictly, Y_b should lie inside the half-plane Y_b > 0 as shown, but half-plane heat conduction solutions for the bed must be evaluated on Y_b, and hence continued outside their domain if Y_b > o. The smoothness conditions on f(x) required for the expansion procedure allow such extensions, so the bed line y = o can be set anywhere in the vicinity of the bed surface. However, changing the bed line y = o for a given profile shape changes f_m and the scale length λ, so that it must be confirmed that this scale change leaves the physical solution invariant. The change is equivalent to a uniform shift of the profile in the y-direction (either direction), so consider

Thus

and if w(x) is any physical variable represented by

on the two scales, then

It will be seen that the subsequent solution is invariant under the transformations shown in Equations (11) and (13), and, by Equation (9), any order of magnitude statement concerning λ applies also to λ_e

An alternative coordinate transformation

proposed by Reference NyeNye (1970) to make the bed surface Y = o requires

so continuing beyond the leading terms gives differential equations of a more complicated form and the biharmonic feature of the theory is lost. Expansion about Y = o of solutions of the original equations is therefore the more appropriate approach.

It is assumed that the regelation process takes place everywhere on the surface producing a continuous thin water layer. The possible shear stress in such a water layer is negligible compared with τ so that one boundary condition for the ice at the bed surface is

where (s, n) denotes the local tangent and normal (inward to the ice) coordinates at a point on the bed surface, with the s-tangent direction inclined at angle θ to the bed line y = o. The condition represented by Equation (16) is an idealization which requires that the water layer is nowhere "pinched-out" or drained away, and which assumes that the ice surface remains in contact with the surface layer; that is, no appreciable cavitation occurs. in terms of the stress tensor σ,

where

The bed surface is everywhere at the pressure-melting point. Let p₀ be a mean pressure level for the bed, and (p₀ T₀ ) a pressure-melting point. Then, for temperatures T close to T₀ there is a linear relation between temperature and pressure-melting point for ice

where C = 0.7 X 10 "⁷ °C m² N^-1 (Reference NyeNye, 1969). Reference KambKamb (1970) points out that the hydrostatic pressure/) is the correct stress in Equation (19). Thus, if T is the temperature in the ice and S the temperature in the bed, temperature continuity at the bed surface gives

if V _n is the normal component of ice velocity on the bed surface (directed into the ice so that a positive value denotes local freezing and a negative value indicates melting), then there is a surface distribution of heat sources with density per unit length of bed LV_n, unit length normal to the plane, where L is the latent heat (Reference NyeNye, 1969) and L = 2.8 × 10⁸ J m^-3. Let N = n/λ be the normal coordinate in (X, Y) coordinates, then the ice and bed temperature fields satisfy the surface flux condition

where k₁, k_b, are the thermal conductivities of the ice and bed respectively. k₁ = 2.0 J m^-1 s^-1 °C^-1 and k^-b = rk_i where r has a range from 1 to 2 for typical bed rocks; the value r == 1.6 appropriate to granite is used in a later calculation.

The stress in the ice is given in terms of the two velocity components V_z V_y, by the constitutive law, so Equations (16), (20) and (21) are four bed surface conditions for V_x, V_y, T, S. On the glacier surface

where p_a is atmospheric pressure. Strictly, the free-surface condition (Equation (22)) is compatible with the restriction y = h only if V_y = O there, but it will be shown that this is valid to any order in e under the weak restriction

since the solution gives exponential decay in Y for V_y. The flow is represented as the sum of a laminar flow satisfying Equation (22) exactly, but not Equation (16), and a perturbation which is not assumed to be small near the bed, only at y = h, so that Equation (22) is satisfied to the required order. The perturbation is constructed as a half-plane solution and Equation (23) is the requirement. Similarly, the ice and bed are treated as half-planes for the heat conduction solutions satisfying Equation (20) and (21) on the bed surface. A geothermal heat flux Q normal to the bed is included by the conditions

A typical value for Q of 4 × 10^-2 J m^-2 s^-1 (Reference PatersonPaterson, 1969) is used Tor later estimates.

3.Flow Equations

The ice is assumed to be an incompressible Newtonian fluid of high viscosity μ. A value for μ of 3 × 10¹² N m^-2 s (Reference NyeNye, 1969) is used in the calculations. Thus

Momentum balance for slow viscous plane flow under gravity requires that

Let

where

and also

where

in these equations u and v are dimensionless velocity components (with unit U_s) while P is a dimensionless pressure (with unit 2μU_s /λ) superposed on a laminar flow which balances the body force (gravity) and which satisfies the upper surface condition (Equation (22)), p₀ is the overburden pressure on a bed line γ = 0 in the laminar flow.

and are dimensionless temperatures (with unit 2μCU _s/λ) in the ice and bed superposed on a uniform flux field satisfying Equation (24). They satisfy the steady heat conduction equations in their respective domains

Here the ice motion is neglected (Reference KambKamb, 1970) and ▽² is the two-dimensional Laplacian in (X, Y) coordinates. Thus, absorbing constants into T₀:

and Equation (20) can be written as the two relations

where

and

D= o if Q = o or if the two conductivities are equal, that is, r = k_b/k_i, = 1. With the values

given earlier, Q/(pgk_iC) = 31.66, and to a good approximation the cos α term in B can be neglected and

This reduction is not appropriate if Q is appreciably smaller than the adopted value, when the full expressions given by Equation (36) must be used. The flux condition, Equation (21), becomes

where

for the adopted values. A value for λ_* of 0.077 m is equivalent to the Reference NyeNye (1969) value for k_* ^-l (in his notation). By Equation (27)

The flow formulae (Equation (25)-(2g)) reduce to

which are the equations for slow viscous plane flow for a velocity field (u, v) and pressure P in the absence of body force and in the ease where the viscosity is 0.5 (Reference LangloisLanglois, 1964). The dimensionless stress Σ with unit 2μU_s /λ is given by

Thus there exist analytic functions φ(ζ),φ(ζ) of the complex variable z = X+iY in h > Y > Y_b such that

Similarly, from Equation (32), there exist analytic functions χι(z), Xb(z) with

4.Flow Solution

Now let u, v,

extend smoothly in Y > o onto Y = o (if Y_b > o anywhere) and to infinity, and let extend smoothly in Y < o onto Y = 0 and to infinity, so that u, v, , S and the stress components are half-plane fields given by the Equation (43)-(47)· Define boundary values Y=o such that

Note that U, V, Θ, Ω are not values of the field variables on the bed surface, but are introduced to allow simple Cauchy integral representations of φ, ψ, χι, χb. it is assumed that they vanish or have sinusoidal behaviour as X → ±∞. Also, we prescribe infinity conditions

and

to be consistent with Equation (33) and the requirement in Equation (22). Then, Equation (43) and (47), subject to Equation (48) and (49), give (Reference MuskhelishviliMuskhelishvili, 1954) :

Once U, V, Θ, Ω are determined, Equation (50)—(53) and (43)-(47) give the physical field variables.

The following results are used repeatedly (Reference MuskhelishviliMuskhelishvili, 1954). If

and W(t) is differentiable and bounded at infinity, then

with extension to higher derivatives, and

if W(t) is continuous and vanishes or behaves sinusoidally at infinity,

where ƒ denotes the Cauchy principal value. Also

where H[W] is the Hubert transform of W(t) (Reference ErdélyiErdélyi, 1954), with the inversion theorem

Note the particular results

which are the only evaluations required when F(X) is a truncated Fourier series, and which show the exponential decay in Y. The first expression of Equation (60) does not have the required conditions at infinity to satisfy the inversion theorem (Equation (59)).

The bed surface conditions (Equation (16), (34), (35), and (38)) can now be expressed in terms of Equation (43)-(47) and (50)-(53); that is, in terms of U, V, Θ and Ω. However, assuming sufficient smoothness, the value of each quantity on Y = _b,() can first be approximated by a truncated Taylor series in Y about Y = 0 to the required order in ε at each X, and then the formulae used to evaluate quantities on Y = 0. If we adopt this process, and omit the lengthy but straightforward algebra using Equation (54)-(58), the bed surface conditions become

The argument X is omitted throughout. if we anticipate the conclusions that U, V, Θ, and Ω are O(ε) or smaller, and that κ(λ/h) is O(ε²) or smaller, then Equation (62)· (65) contain explicitly the terms required to derive expansions of U, V, Θ and Ω to O(ε³); that is, an approximation neglecting terms O( ε⁴) compared with a leading term O(ε). Thus, for a smooth bed, the error may be only 1% for ε = O.2. Here terms to O(ε ²) only are derived explicitly. it is supposed that B and D are order unity in the balance, but later calculation with the adopted physical parameters shows that these geothermal flux terms contribute little to the O(ε) terms.

Consider the balance for λ = ?(λ_*), that is

Alternative balances for smaller and larger ω have been obtained explicitly, and have been shown to be appropriate limits of the solution based on Equation (66). Let

with similar expressions for V, Θ and Ω, and compare coefficients of ε⁰, ε¹, ε², in Equation (62)-(65). The leading term of Equation (62) gives

and by Equation (59) and (6o),

since a non-zero constant U₀ leads to finite φ{?) and ψ(?), as Y, → ∞, which violates Equation (49). By Equation (64), (65) and (63) in turn,

and eliminating (Θ)₀using Equation (58) and (59),

Both solutions of Equation (71) are unbounded at one limit, X → ± ∞, so

With Equation (69) and (72), the ε term of Equation (62) is analogous to Equation (68), giving

Now Equation (64), (63) and (65) give

where

and

for the adopted values (r = 1.6). λ_* = 0.088 m is the Reference KambKamb (1970) definition of the critical length scale. Both complementary functions of Equation (75) are unbounded at one limit, X > → ± ∞, and V₁ , is given by the bounded particular integral, simply a repeated quadrature once F' and H[F'] are known.

Now the ε² term of Equation (62) gives

H[U₂ ] and H[V₁”] are periodic, but the product FH[V₁”] is in general the sum of a constant term Γ and a periodic term W(X), so the solution of Equation (78) is

This is demonstrated clearly by the subsequent sinusoidal bed calculation. Thus,

By Equation (75) Γ contains a term (1 — κ), so an explicit relationship between κ, λ/h and ε² is obtained for a given F(X), and in turn we obtain expressions for U_s and U_b/U_s where by Equation (27) and (69):

The basal sliding velocity U_b is defined as the term of order unity in V_x on Y = 0. It will be shown in Figure 3 that the terms in B and D do not contribute to Γ, and hence do not affect U_s and U_b, and furthermore that 0 < κ<1 as expected. It should be noted that the second-order velocity term U₂ and the basal sliding velocity (of order unity) are determined together in an O (ε²) balance.

Finally, to complete the second-order expansions, using Equation (64), (63) and (65),

The latter has a particular integral for V₂ , which is bounded and periodic. The product terms in Equation (82) allow constants to appear in Θ₂ and Ω₂. These represent second-order corrections to T₀ in the values of T and S on Y - 0. Recall that T₀ was defined as the melting temperature corresponding to pressure p_u

5.Sinusoidal Bed Profile

The solution is now given explicitly for

with F(X) = cos X covered by the shift X-> X + π/2. The contributions of the harmonics in a truncated Fourier series will be indicated. Only the Hubert transforms (Equation (60)) are required. We set

as a result of which, Equation (75) becomes

with the bounded particular integral

Hence

so that

and

Thus κ, and hence U_b and U_s , are independent of A for F(X) = sin X. This applies for any combination a sin X + b cos X which becomes a sine function with a phase shift. If further harmonics are added, the cross products between different harmonics arising in FH[V₁ ”] contribute no constant terms, so A cannot enter into Equation (89). The restriction of Equation (23) is met provided that κ > 0(ε²) for ≥ 0(1), thus implying an upper bound to ? (see Equation (91)).

From Equation (28) and (81) and the leading term of Equation (89},

For a fixed ratio of h/ε², the minimum value of U_b is given by:

and this occurs when λ = λ_*. So, for λ smaller or larger than λ_*, the basal sliding velocity increases. Figure 2 shows the variation of U_b (in SI-units) with log (λ/λ_*)) for the physical data given earlier. At fixed values of λ, U_b, increases with h/ε². Figure 3 shows the variation of U_b /U_s with log (λ/λ_*) for different values of hε². For example, hε² = 4 applies for h = 100 m and € = 0.2, while hε² = 0.1 applies for h = 10 m and ε = 0.1. However, the Newtonian approximation could have a serious effect on any quantitative predictions of U_b and U_s, whereas the main purpose was to show that a flow field, which determines U_b and which is compatible with the gravity drive, pressure-melting, and regelation boundary conditions, exists. Alternatively, given U_b , U_s and h, Equation (91) and (92) determine the bed parameters ε and λ.

Fig. 2. Variation of basal sliding velocity U_b, with bed. Scale λ.(SI-units).

Fig. 3. Variation of the ratio of basal sliding velocity to surface velocity with bed scale λ for different values of hε² (SI-units).

Substitution of the above expressions indicates that the right-hand side of Equation (83) is zero, and hence the bounded solution is

The temperature terms are now given by Equation (74) and (82). The velocity components of the bed line are, to second order,

Using the estimates contained in Equation (37), (39), (85) and (89),

Taking the extreme values h = 10 m, α = 0.04 and ε = 0.2, the maximum value of A is 0.05, so thai the contribution made by the geothermal flux is small. We may recall that B and hence A, are even smaller if Q = 0, by Equation (36).

The substitution of Equation (95) into Equation (50) and (51) and the use of Equation (61) gives, to second order,

which determine the velocity and stress fields through Equation (43) (46) and Equation (27)-(30), and which show the exponential decay in Y. In particular

We may note that

and hence the pressure field becomes

when A is neglected.

On the bed surface

and since the present solution is valid only when the ice boundary coincides with the bed surface (that is, no cavitation occurs) we require p_b ≥ 0 which implies that

The factor (1 + p_a/(pgh cos α)) is approximately equal to one for h ⩾ 100 m, and even for a thin glacier h = 10 m (Equation (102)) only allows α Wlsims; ε. So, cavitation is predicted for moderate bed inclination. Reference KambKamb (1970) infers an instability criterion of α gt; ε from an argument involving the value of the absolute bed slope, but this appears to be invalid unless the resisting "pressure drag" is eliminated from the lee flanks of the undulations by "total cavitation" there. While Equation (102) shows that cavitation does occur for α gt; ε (for a sinusoidal bed on the Newtonian approximation), it is still reasonable to expect that pressure over the contact section of lee flanks will be sufficient to provide the necessary drag, but a full solution which incorporates cavitation is needed in order to confirm this.

6.Concluding Remarks

We now see that slow viscous plane flow over an inclined wavy bed has a bounded solution consistent, at least to the order of ε², with the perturbation on the laminar flow decaying rapidly with height. This is in contrast to calculations of half-plane flow which are subjected to the classical zero velocity condition at the bed-line. The inclusion of gravity in the momentum balance provides the driving mechanism and also defines the mean drag directly. Thus, the unbounded terms arising from the artificial driving shear stress of Nye do not arise. The terms in “ agree with the first-order expansions of Nye and Kamb, but, whereas Nye derives a relationship between the mean drag and the basal sliding velocity from the first-order solution, the present solution shows how the surface and basal velocities are governed by the second-order balance and are given directly. it is shown that a uniform geothermal flux does not influence the basal-sliding velocity, and makes only a small contribution to the flow field.

The application of boundary conditions to the bed surface instead of the half-plane boundary allows direct expansions beyond the first power in ε The balance equations have been given up to terms in ε³ so that a solution for moderate slopes (ε ~ 0.2) could be obtained. In the dimensionless expansion analysis, coefficients should remain of order unity if the higher derivatives of the bed profile are also small. However, the complex variable method of conformai mapping offers a direct approach to the problem of flow over a hump of finite slope.

Acknowledgement

This investigation was pursued in connection with a Natural Environment Research Council Grant GR3/2680 "Flow of glaciers over deformable materials" which i hold jointly with Dr G. S. Boulton of the School of Environmental Sciences, University of East Anglia.

MS. received 20 September 1975 and in revised form 8 January 1976