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 pa 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 Ub, and surface velocity Us, in terms of h, α, and the bed geometry for given ice and bed properties. An illustration for a sinusoidal bed profile shows that Ub is very sensitive to the wavelength of the undulation. The variations of Ub 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 ε2 obtained explicitly, though the expansions given determine valid solutions to order ε3 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 ε2 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 ε2 terms, so it is essential to construct a consistent expansion scheme to order ε2. 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, Us is the surface velocity at y = h, and (Ub 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,
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 fm 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 = Yb, to a required order in ε. Strictly, Yb should lie inside the half-plane Yb > 0 as shown, but half-plane heat conduction solutions for the bed must be evaluated on Yb, and hence continued outside their domain if Yb > 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 fm 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 p0 be a mean pressure level for the bed, and (p0 T0 ) a pressure-melting point. Then, for temperatures T close to T0 there is a linear relation between temperature and pressure-melting point for ice
where C = 0.7 X 10 "7 °C m2 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 LVn, unit length normal to the plane, where L is the latent heat (Reference NyeNye, 1969) and L = 2.8 × 108 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 k1, kb, are the thermal conductivities of the ice and bed respectively. k1 = 2.0 J m-1 s-1 °C-1 and k-b = rki 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 Vz Vy, by the constitutive law, so Equations (16), (20) and (21) are four bed surface conditions for Vx, Vy, T, S. On the glacier surface
where pa is atmospheric pressure. Strictly, the free-surface condition (Equation (22)) is compatible with the restriction y = h only if Vy = 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 Vy. 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 × 1012 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 Us) while P is a dimensionless pressure (with unit 2μUs /λ) superposed on a laminar flow which balances the body force (gravity) and which satisfies the upper surface condition (Equation (22)), p0 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 ▽2 is the two-dimensional Laplacian in (X, Y) coordinates. Thus, absorbing constants into T0:
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 = kb/ki, = 1. With the values
given earlier, Q/(pgkiC) = 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μUs /λ is given by
Thus there exist analytic functions φ(ζ),φ(ζ) of the complex variable z = X+iY in h > Y > Yb 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 Yb > 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(ε2) or smaller, then Equation (62)· (65) contain explicitly the terms required to derive expansions of U, V, Θ and Ω to O(ε3); that is, an approximation neglecting terms O( ε4) compared with a leading term O(ε). Thus, for a smooth bed, the error may be only 1% for ε = O.2. Here terms to O(ε 2) 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 ε0, ε1, ε2, in Equation (62)-(65). The leading term of Equation (62) gives
and by Equation (59) and (6o),
since a non-zero constant U0 leads to finite φ{?) and ψ(?), as Y, → ∞, which violates Equation (49). By Equation (64), (65) and (63) in turn,
and eliminating (Θ)0using 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 V1 , is given by the bounded particular integral, simply a repeated quadrature once F' and H[F'] are known.
Now the ε2 term of Equation (62) gives
H[U2 ] and H[V1”] are periodic, but the product FH[V1”] 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 ε2 is obtained for a given F(X), and in turn we obtain expressions for Us and Ub/Us where by Equation (27) and (69):
The basal sliding velocity Ub is defined as the term of order unity in Vx 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 Us and Ub, and furthermore that 0 < κ<1 as expected. It should be noted that the second-order velocity term U2 and the basal sliding velocity (of order unity) are determined together in an O (ε2) balance.
Finally, to complete the second-order expansions, using Equation (64), (63) and (65),
The latter has a particular integral for V2 , which is bounded and periodic. The product terms in Equation (82) allow constants to appear in Θ2 and Ω2. These represent second-order corrections to T0 in the values of T and S on Y - 0. Recall that T0 was defined as the melting temperature corresponding to pressure pu
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 Ub and Us , 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[V1 ”] contribute no constant terms, so A cannot enter into Equation (89). The restriction of Equation (23) is met provided that κ > 0(ε2) 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/ε2, the minimum value of Ub is given by:
and this occurs when λ = λ*. So, for λ smaller or larger than λ*, the basal sliding velocity increases. Figure 2 shows the variation of Ub (in SI-units) with log (λ/λ*)) for the physical data given earlier. At fixed values of λ, Ub, increases with h/ε2. Figure 3 shows the variation of Ub /Us with log (λ/λ*) for different values of hε2. For example, hε2 = 4 applies for h = 100 m and € = 0.2, while hε2 = 0.1 applies for h = 10 m and ε = 0.1. However, the Newtonian approximation could have a serious effect on any quantitative predictions of Ub and Us, whereas the main purpose was to show that a flow field, which determines Ub and which is compatible with the gravity drive, pressure-melting, and regelation boundary conditions, exists. Alternatively, given Ub , Us and h, Equation (91) and (92) determine the bed parameters ε and λ.
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 pb ≥ 0 which implies that
The factor (1 + pa/(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 ε2, 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 ε3 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