Introduction
The mode of flow of water at the beds of temperate glaciers has been of considerable interest to glaciologists, particularly because of apparent relationships between subglacial hydrology and temporal variations in sliding velocity (e.g.Reference Weertman Weertman, 1972, p. 288–89; Reference Hodge.Hodge, 1974,Reference Hodge1979; Reference Röthliskberger, Röthlisberger, Iken and SpringRöthlisberger and others, 1979). A number of recent field studies of existing glaciers have yielded important data on subglacial hydrology (Reference MathewsMathews, 1964,Reference Mathews 1973; Reference MeierMeier, 1965; Reference Stenborg Stenborg, 1969, Reference Stenborg1973; Reference HodgeHodge, 1976, Reference Hodge1979; Reference Engelhardt Engelhardt, 1978; Reference EngelhardtEngelhardt and others, 1978; Reference Iken, Iken, Flotron, Haeberli and RöthlisbergerIken and others, 1979; Reference Röthliskberger, Röthlisberger, Iken and Spring Röthlisberger and others, 1979). Water-level soundings in bore holes drilled to glacier beds clearly indicate the existence of channels in which the water pressure may be substantially less than the overburden pressure. Fluctuations in bore-hole water levels have been attributed to variations in both melt-water input and the geometry of the subglacial drainage system. Furthermore, the probability of bore holes connecting to subglacial water channels may vary widely from year to year at any point on a glacier. These studies have also provided estimates of the amount of water stored within glaciers and the speed at which water moves through the ice.
Studies of recently deglaciated rock surfaces have also provided important data on the geometry of subglacial drainage networks. Reference Walder and HalletWalder and Hallet (1979)and Reference Hallet and Anderson.Hallet and Anderson (in press) mapped systems of interconnecting channels and cavities comprising about 20–40% of the former beds of two small cirque glaciers. Variation in the subglacial hydraulic network through time was clearly indicated by the fact that many channels incised into the rock by turbulent subglacial melt water were also coated with CaCO3 precipitate that forms only where the rock is in intimate contact with sliding ice (Reference HalletHallet, 1979). The mapping by Reference Walder and HalletWalder and Hallet (1979) further revealed 5–10 m wide zones, paralleling the former direction of ice flow, which were thoroughly striated and nearly devoid of the precipitate. Apparently, the bed in these zones had been covered by a water sheet accommodating a water flux sufficient to inhibit the build-up of the solute concentration to the saturation point.Reference Hallet Hallet (1979) has estimated from the size distribution of rock fragments in the subglacial precipitate, that over most of the bed of the glacier studied, the water sheet was typically no thicker than several tens of micrometers, possibly thickening at times to c. 100 μm.
The picture of the subglacial drainage system arising from these various studies is similar to that envisioned by Reference NyeNye (1973): a thin water sheet over most of the bed, co-existing with but relatively independent of discrete channels that drain most of the melt water. Furthermore, both channelized and sheet flow are variable in time and space. Theoretical models should account for these observations. Significant progress toward the development of such models has been made by Reference Röthlisberger.Röthlisberger (1972) and Reference ShreveShreve (1972), who showed independently that channelized flow tends to become concentrated, with larger channels growing at the expense of smaller ones, forming an arborescent network. Both of these authors argued that the stable form of subglacial water passages should be tunnel-like rather than sheet-like. However, Reference WeertmanWeertman (1972) showed that discrete channels at the glacier bed might not be efficient melt-water collectors, due to the effective linearization of the creep law for ice in the particular case of tunnel closure when the relatively high shear stress at the bed is explicitly considered. Weertman showed that this linearization of the creep law could result in pressure gradients that would actually drive water away from channels at the glacier bed. He therefore concluded that primary drainage of melt water at the beds of temperate glaciers must be through a widespread, relatively thick sheet of water.Reference Nye Nye’s (1976) study of a jökulhlaup, however, has led him to conclude otherwise. Although Nye did not explicitly determine the temperature distribution, and hence the heat flow, in a sheet of water at the bed of a temperate glacier, he argued that heterogeneous heat production in a water sheet of variable thickness should lead to flow localization. Nye’s argument serves as a starting point for the analysis presented here.
In this paper, I have analyzed in detail the stability of sheet flow of water at the bed of a temperate glacier with respect to perturbations in the sheet thickness. I also discuss the implications of the results for the well-known water-lubrication theory of glacier surges (Reference WeertmanWeertman, 1962, Reference Weertman1964,Reference Weertman 1966, Reference Weertman.1969).
2. Outline of the analysis
Sheet flow of water beneath a temperate glacier tends to be unstable with respect to variations in the sheet thickness, i.e. such variations tend to become enhanced. Physically this effect arises from the greater viscous heat dissipation due to water flow, and concomitant higher melting rate, where the sheet is thickest, as compared to where it is thinnest. In idealized sheet flow, two processes tend to counteract the instability: cross-stream heat flow within the sheet and “plastic” sagging of the ice. Furthermore, incipient channels may be destroyed as the glacier slides over protuberances on its bed (Reference NyeNye, 1973). The central issue to be addressed is therefore the rate of growth of perturbations in the water sheet.
In section 3, I present solutions for the fluid flow and heat transfer within a subglacial water sheet that varies slightly in thickness along a direction normal to the flow. The analysis is rigorous to first order in the small thickness perturbation. Adoption of the sheet-flow approximation is tantamount to assuming that the glacier bed is smooth on a scale large compared with the sheet thickness. Formal incorporation of the effects of glacier bed irregularities and subglacial drift would greatly complicate the model, but would be unlikely to alter our picture of the fundamental physics involved.
Sagging of the ice is analyzed in section 4 by considering the ice-water contact to be the interface between two linear viscous fluids and studying the relaxation behaviour of perturbations on an otherwise planar interface, following the method of Reference FletcherFletcher (1977). In section 5, I discuss how the roughness of real glacier beds may be incorporated into the model in an approximate fashion, as well as the implications of my results forReference Weertman Weertman’s (1962, Reference Weertman1964, Reference Weertman1966, Reference Weertman.1969 water-lubrication theory of glacier surges.
3. Fluid Flow and Heat Transfer in the Water Sheet
In this section, I derive expressions for the flow velocity and temperature distribution in a subglacial water sheet of variable thickness. Several simplifying assumptions are made about the geometry. Water flow is assumed to be steady, one-dimensional, and parallel to the ice flow. The bed is considered to be essentially planar, with the surface slope α of the ice providing a constant pressure gradient driving the flow. The water pressure is assumed to be equal to the ice overburden pressure p. In reality, sliding of the ice over bedrock protuberances will result in regions of high pressure where the water sheet is largely squeezed out. Water will then tend to flow around these regions, resulting in a slightly lower average pressure gradient in the sheet (Reference WeertmanWeertman, 1972, p. 311–12), as well as two-dimensional flow. Again, this effect is unlikely to alter the results of the simple model fundamentally.
The coordinate system and idealized water-sheet geometry are depicted in Figure 1, with x representing distance along the bed in the flow direction, y perpendicular to x and along the bed (the “lateral” direction), and z perpendicular to the bed. The glacier bed, assumed impermeable, is at z = 0. The water-ice interface is described by the expression
where h is the average thickness of the sheet and ε ⪡1. The amplitude Ȧ of the perturbation therefore equals εh. The Navier- Stoke equations for steady flow reduce to (Reference Bird, Bird, Stewart and LightfootBird and others, [c 1960], p. 80)
where u is the flow velocity, P g = –∂p /∂x is the pressure gradient driving the flow, the sign chosen as to make P g > 0, and ηw is the viscosity of water.
Fig. 1. Coordinate system and idealized water-sheet geometry for the theoretical analysis.
I now assume that the flow velocity can be expressed as a mean flow plus a small perturbation:
Substituting Equation (3) into Equation (2) and separating terms of O(1) and O(ε) leads to the two differential equations
The boundary conditions on the flow are
These are to be satisfied to O(ε). The following velocity distribution in the sheet, which satisfies these boundary conditions, is derived in Appendix B:
This expression is used in computing the temperature field in the water sheet. The thermal energy equation for steady-state, incompressible, one-dimensional flow is (Reference Bird, Bird, Stewart and LightfootBird and others, [c 1960], p. 315)
where T is the water temperature; ∇2 is the three-dimensional Laplacian operator in Cartesian coordinates, p w, c w, and kw are the density, constant-volume specific heat, and thermal conductivity, respectively, for water at the melting point.
In Appendix C, it is shown by dimensional analysis that both down-stream and lateral heat conduction, corresponding respectively to the terms k w ∂ 2 T/∂x 2 and k w ∂ 2 T/∂y 2, are negligible in the overall heat balance in nearly all instances, including the interesting case of a sheet thick enough substantially to affect bed roughness. However, the small, but finite, amount of lateral heat conduction might be expected to affect the growth rate of perturbations; hence the term k w ∂2 T/∂y 2 will be retained in the thermal energy balance.
One further simplification is made to permit an analytic solution to Equation (6). The temperature of the water in a thin sheet is expected to differ only very slightly from the pressure-melting temperature, hence the approximations
where T m is the pressure-melting temperature and c t = –∂T m/∂p, the choice of sign making c t > 0. Since a constant pressure gradient has been assumed, ∂T/∂x is constant in this approximation. Substituting Equations (8) into Equation (7) and neglecting down-stream conduction, the thermal energy equation reduces to
where γ= p w c w c t = 0.316 is a dimensionless constant that arises from the pressure-melting behavior of ice (Reference Röthlisberger.Röthlisberger, 1972).
The temperature field is now decomposed into a mean and a small perturbation:
The boundary conditions on the temperature field are
where qG is the geothermal heat flux at the bed. Equation (9) can be transformed into two equations to be solved independently for T0 and T1, subject to the boundary conditions of Equations (11). Details of the rather lengthy analysis are given in Appendix D. The quantity of greatest interest for the stability analysis is the heat flux from the water sheet into the basal ice, –k ∂T/∂n where n is the local upward normal from the ice-water interface. Applying the chain rule for differentiation and noting from Equation (1) that the slope ∂z0/∂y of the interface is O(ε), it follows that, to O(ε),
The heat flux into the basal ice, denoted QZ0 , is given by Appendix D:
where I have introduced the approximation kh⪡ 1, i.e. the wavelength (= 2π/k) of the perturbation is large compared with the sheet thickness. For such values of kh, the functional form of the heat flux is most readily understood.
In Equation (13), the heat flux has been separated into three bracketed members. The first is the mean heat flux, composed of terms representing geothermal heat and viscous dissipation. The factor (l–y)≈2/3 corrects for the net down-stream advection of heat that maintains the water at the pressure-melting temperature.
The second bracketed member of Equation (13) is the locally enhanced viscous heat production in the water sheet, again corrected for down-stream advection. Heat flow into the ice, and therefore melting, is greater in the thick parts of the sheet (sin ky > 0) than in the thin parts (sin ky <0). This is now clearly seen to be the source of instability.
Finally, the third member of Equation (13) accounts for heat conducted laterally from thick to thin parts of the sheet, as a result of the warping of the isotherms in the sheet, relative to the unperturbed case. This quantity is seen to be O(k 2 h 2) and thus of negligible importance for kh ⪡ 1, in agreement with the scaling arguments of Appendix B.
The rate at which perturbations grow due to melting can now be examined. The melting rate of the basal ice is Qz0/p i L, where p i is the density of ice and L is the heat of fusion. Neglecting the mean heat flux, which causes equal amounts of melting everywhere on the ice-water interface and does not affect perturbations, the speed w m at which the interface moves due to melting is found to be
in the next section, the speed at which the ice-water interface moves due to sagging of the ice is determined. The sum of that speed and w m will then determine the growth rate of perturbations.
4. Sagging of the Ice and Growth Rate of Perturbations
The relaxation of the perturbed ice-water interface can be analyzed by treating the interface as one between two fluids of different densities and viscosities. This approach is widely used in the theory of folding of layered materials; I apply here the particular method developed by Reference FletcherFletcher (1977). In this formulation, the fluids are assumed incompressible and quasi-static, i.e. accelerations are ignored. The rheology is taken as linear viscosity. The choice of rheology is particularly suitable for this model, because for incipient channels at the glacier bed the difference between ice and water pressures will not be large compared to the shear stress, thus effectively linearizing the creep law (Reference WeertmanWeertman, 1972, p. 299–306). Stresses and velocities are then each decomposed into a mean and a small perturbation. In this particular model, the mean stress is hydrostatic pressure, while the mean velocities are identically zero. Velocity components v and w are associated with the relaxation process and correspond to the y- and z-directions. respectively. Boundary conditions on stresses and velocities are satisfied to first order in the slope of the perturbed interface, i.e. to O(kA). where kA ⪡1 is assumed.
In practice, the effective viscosity of ice is so much larger than that of water, that the analysis reduces to a determination of the speed at which perturbations relax on a free surface. The quantity of interest is w p, the speed at which the interface moves due to viscous sagging. From Appendix E, Equation E-3,
where g is the acceleration due to gravity and η i is the effective viscosity of ice.
The mean rate of melting of the basal ice maintains an average sheet thickness of h; heterogeneities in melting and sagging rates cause the amplitude of the thickness perturbation to change quasi-statically, but without any change in the mean thickness of the water sheet. Adopting this view, points on the interface z = z 0 can be considered as remaining on the interface. Mathematically, this is expressed as (Reference FletcherFletcher, 1977, p. 600)
where D/Dt = (∂/∂t) + v(∂/∂y) + w(∂/∂z) is the total time derivative. Expanding Equation (16) yields
Substituting Equations (1), (14), and (15), recalling thatA = hε, and noting from Appendix E that v is of O(kA), I find to O(kA), after rearranging:
where
and kh⪡ 1. Table 1 lists values of l/τ1, 1/τ2, and 1/τ3 for several choices of P g , h, and k; the values of the other physical constants are in Appendix A. It is clear that l/τ3 is negligible in all instances and that 1/τ 2 ⪡ 1/τ1 unless A and k are very small. Neglecting these cases and noting the restriction kh ⪡ 1, Equation (18) can be closely approximated as
The amplitude of small perturbations in the water sheet therefore grows exponentially with time constant τ1. The strength of the instability increases with increasing P g and h. Figure 2 shows τ1, as a function of h for three reasonable choices of P g , the term k 2 h 2 being neglected.
Fig. 2. Characteristic growth and decay times (τ1 and τd respectively) for water-sheet perturbations.
Table I. Time constants describing rate of growth of water sheet perturbations
5. Effects of Bed Roughness and Implications for Glacier Surging
The analysis above demonstrates that sheet flow on a planar glacier bed would be unstable with respect to perturbations in sheet thickness. The roughness of actual glacier beds will alter this simple result in two ways, however. First, perturbations in the water sheet, also referred to as “incipient channels”, may be destroyed as the glacier slides over protuberances on the bed (cf.Reference NyeNye, 1973). Second, sub-glacial cavities on the lee sides of bed protuberances may capture significant amounts of melt water and limit the thickness of water sheets (Reference Walder and HalletWalder and Hallet, 1979; Reference Hallet and Anderson.Hallet and Anderson, in press).
The rate at which incipient channels in the basal ice encounter bed protuberances must generally increase as sliding velocity and bed roughness increase. The criterion for sheet stability used herein is that the sheet will be considered quasi-stable when the average time between encounters with bed protuberances, hereinafter referred to as the “decay time” τd, exceeds the characteristic time for growth of sheet perturbations.
A rigorous analysis of the rate at which incipient channels are destroyed by encounters with bed protuberances would involve coupling glacier sliding physics with that of water flow in a “sheet” of complex geometry and spatially variable pressure gradient. Such is beyond the scope of this paper. Nonetheless, a reasonable first approximation to the solution of this problem may be found by adopting a glacier bed model of small hemispheres distributed on a plane. This model was used by Reference LliboutryLliboutry (1978) in his formulation of glacier sliding theory. As the glacier slides, incipient channels in the basal ice will encounter these bed irregularities and may be destroyed if the irregularities are neither very small nor very large compared to R⋆ ., the “transition obstacle size” in Reference LliboutryLliboutry’s (1978) glacier sliding theory. Very small irregularities would be submerged by the water sheet, whereas the glacier would slide over large protuberances predominantly by plastic deformation with very little regelation. In the latter case, flow lines in the ice would be very nearly parallel to the bed; hence, the form of the incipient channel could be preserved.
Following Reference LliboutryLliboutry (1978, p. 152), I assume a “non-dimensional” bed, such that the fraction of the bed covered by hemispherical bed irregularities with radii in the range (R, R + dR) is equal to μ dR/R where μ, a constant, is a measure of the bed roughness. The fraction of the bed covered by all irregularities is thus μ ln (R max/R min), where R max and R min are, respectively, the maximum and minimum radii of bed irregularities. This fraction is assumed to be ⪡ 1. Lliboutry suggests R min = 1 μm, Rmax = 10 m, hence In (Rmax/Rmin)≈16 and μ must be ⪡ 1/16. Due to the logarithmic dependence, R max and R min may vary significantly from these suggested values without appreciably affecting the constraint on μ.
As discussed above, bed irregularities with radii near R ⋆ are most “effective” at destroying incipient channels. I will use the following form for the “relative efficiency”:
which is quite similar to the expression found by Reference WattsWatts (unpublished) for the drag on a sphere moving through temperate ice. The chosen functional form of E has the properties that E(R⋆) = 1, and that E approaches zero as R approaches zero or infinity.
The number of bed irregularities per unit area of the bed, with radii in the range (R, R + dR), is (Reference LliboutryLliboutry, 1978, p. 152)
hence, the “effective” number v e of bed irregularities per unit bed area, i.e. the number that may actually destroy an incipient channel, is
where the lower limit of integration reflects the size of the smallest bed irregularities that “block” an incipient channel of width W.
The average number of the “effective” irregularities along the flow direction in a rectangular bed area of length l along the flow direction, width W, is velW. Because the average spacing of these irregularities is l/v e lW, the average time between encounters with such irregularities can be expressed as
where I have evaluated the integral in Equation (22). Reference Lliboutry Lliboutry (1978, p. 152) suggested that R⋆ = 0.16U −1/2, for R⋆ in meters, U in meters per year. Although the coefficient of 0.16 is probably incorrect, due to the fact thatReference Lliboutry Lliboutry’s (1978) glacier sliding model does not fully satisfy thermal boundary conditions at the glacier bed (recognized by Reference LliboutryLliboutry (1979, p. 80)), increasing or decreasing this coefficient by a factor of two can be shown not to affect the results presented herein substantially. Hence, using Lliboutry’s expression, Equation (23) becomes
The width W of an incipient channel can be identified with the wavenumber k of the water sheet perturbation by the relationship W= πk−1 ; i.e. W is the perturbation half-wavelength. The analysis for the growth of perturbations is valid for kh⪡ 1; hence, the width of incipient channels that can be examined by the present analysis is constrained by the relationship W⪢πh. A reasonable minimum value for W is therefore 30h. Using this value, I can rewrite Equation (24) as
where h is expressed in meters, U in meters per year, τd in years, and μ = 0.01 in accord with the constraint discussed above. This relationship has been plotted in Figure 2 for three choices of U. The values of τd computed from Equation (25) are minimum values, because τd increases with W. The criterion for water-sheet stability used here is that τd < τ1 for the specified values of Pg and U; hence, by using Equation (25) for τd, the probability of sheet stability is maximized.
It appears from Figure 2 that, for reasonable choices of Pg and U, the stability criterion can be met only for h <c. 1–4 mm. Sheets of greater thickness will tend to be unstable.Reference WeertmanWeertman (1962, Reference Weertman1964,Reference Weertman1966, Reference Weertman.1969) argued that a water sheet thick enough substantially to reduce effective bed roughness will cause an increase in sliding velocity of one to two orders of magnitude, i.e. a glacier surge. He approximated bed roughness by cubes, with those of side Λ being most important in causing drag on the sliding ice. The homologous quantity in (Reference LliboutryLliboutry’s 1978) theory is R ⋆, which can be expressed in terms of U. If a sheet of thickness h ≥ R ⋆ is stable to thickness perturbations, then Weertman’s mechanism might operate.
Table II gives values of h max/R ⋆, where h max is the maximum stable thickness, for several choices of P g and U. In all cases examined, encompassing the great majority of actual glaciers, h max/R ⋆ is much less than unity; the largest values of h max/R ⋆ are found when P g is relatively small and U relatively large, i.e. for gently sloping, fast-moving glaciers. Hence, the conclusion reached from the present model is that the Weertman water-lubrication mechanism is unlikely to be effective for surge initiation. However, this mechanism may be effective at maintaining surges once high sliding velocities are achieved. Caution must be taken with this conclusion, however. If bed irregularities of width R < 1/2W can cause destruction of incipient channels, then τd could be reduced considerably from that given by Equation (25) with the result that thicker sheets could be stable.
Table II. Ratio of the maximum stable sheet thickness hmax to the transition obstacles size R*
An independent and fundamentally important consideration in evaluating the water-lubrication theory of surges is whether sufficient melt water is ever available at the bed to cause thick water sheets. The thickness of such a sheet at the glacier bed may be computed by assuming that uniform melting of thickness M per unit time occurs everywhere on the base of the glacier, and that all melt water flows at the bed. The sheet thickness is then (Reference Weertman.Weertman, 1969, p. 953)
where x = 0 is the point farthest from the glacier terminus. Figure 3 shows h as a function of x for several values of M and P g, the latter also expressed as surface slope α. For typical values of geothermal heat flux (c. 0.05 Jm−2 s−1), a bed shear stress of c. 1 bar, and sliding velocity U up to several tens of m a−1, M = c. 15 mm a−1. Figure 3 shows that with such a melting rate, h ≤ 2 mm for glaciers less than c. 100 km long, regardless of surface slope. Such a value of h is less than estimates of the height of the most important bed-roughness elements. Thus, “normal” melting rates are very unlikely to produce water sheets thick enough to substantially reduce bed roughness.
Fig. 3. Thickness of water sheet as a function of distance along the ice-flow direction, assuming all melt water flows in the sheet.
In order for h to reach even 5 mm, a likely lower limit for the “controlling obstacle size”, M must be considerably larger than 15 mm a−1. Studies of intergranular vein structures in temperate glacier ice (Reference Raymond and HarrisonRaymond and Harrison, 1975), and of ice chemistry (Reference Berner, Berner, Stauffer and OeschgerBerner and others, 1978), suggest that as much as c. 100 mm a−1 of melt water can percolate to the glacier bed through intergranular veins. However, even for an effective melting rate of M = 120 mm a−1, an eight-fold increase from the “normal” melting rate, h will be increased by a factor of only two and will therefore still be less than 5 mm for glaciers of lengths less than c. 100 km.
It appears that other distributed sources of large amounts of melt water are necessary if the subglacial water sheet is to become thick enough to initiate surges. As shown in Figure 3, an effective melting rate of 2 m a−1 could produce a sheet of thickness c. 5 mm for glaciers with lengths c. 10–50 km and surface slopes of c. 0.5–4.0. A number of surging glaciers have these dimensions (Reference Meier and PostMeier and Post, 1969). It is clear that water flows through englacial and subglacial conduits, as well as through intergranular veins (Reference MathewsMathews, 1964; Reference Röthlisberger.Röthlisberger, 1972; Reference ShreveShreve, 1972; Reference Weertman.Weertman, 1972;Reference Nye Nye, 1973, Reference Nye and Frank1976;Reference Fletcher Vivian and Zumstein, 1973; Reference Hodge.Hodge, 1974, Reference Hodge1976, Reference Hodge1979; Reference FletcherRaymond and Harrison, 1975;Reference EngelhardtEngelhardt, 1978;Reference Engelhardt Engelhardt and others, 1978; Reference Röthlisberger.Röthlisberger and others, 1979; Reference Walder and HalletWalder and Hallet, 1979), but estimates of the magnitude of such flows are not available. Furthermore, it is possible that melt water reaching the bed through such conduits would remain channelized, rather than spreading out (Reference ShreveShreve, 1972). The fact that bore holes drilled to glacier beds often encounter channels suggests that such melt water does not join a subglacial water sheet. Considerably more information is needed about the hydrology of temperate glaciers before it can be established that major distributed water sources can actually supply a subglacial water sheet.
A final consideration, related to the above discussion of melt water supply, is the amount of water flow through subglacial cavities.Reference Walder and Hallet Walder and Hallet (1979) andReference Hallet and Anderson. Hallet and Anderson (in press) have shown that cavities may comprise as much as 20–40% of the ice-rock interface beneath two small cirque glaciers. Hallet and Anderson have calculated that the potential water storage in an exhumed cavity network at Castleguard Glacier, Alberta, would be equivalent to a sheet of water 63 mm thick, if all cavities were filled. Reference HodgeHodge (1974) estimated that the water stored at the beds of Nisqually Glacier and South Cascade Glacier, later released during jökulhlaups, could amount to an equivalent water layer approximately 1 m thick. Furthermore, because the subglacial cavities are regions of relatively low water pressure, they tend to act as continual “sinks” for subglacial melt water. The total of the water stored in subglacial cavities, plus the through-flow of water in the cavity-channel network (Reference Walder and HalletWalder and Hallet, 1979), could amount to a significant fraction of the total annual melt-water production. The water sheet thickness shown in Figure 3 must therefore be considered an upper limit, applicable only in the absence of numerous subglacial cavities. Hence, the presence of water-filled subglacial cavities reduces the possibility that the water lubrication mechanism for glacier surges can operate.
6. Conclusions
I have presented a model for the stability of sheet flow of water at the base of a temperate glacier, with respect to perturbations in sheet thickness. Sheet flow is nearly always unstable on planar glacier beds, but the roughness of real beds may result in quasi-stable sheet flow for sheets of thickness no greater than c. 4 mm. However, the choice of model for bed roughness has a strong influence on the values determined for maximum stable sheet thickness.
Quasi-stable subglacial water sheets are apparently not thick enough to cause drastic reductions in bed roughness and thereby initiate surges, but may cause some increase in sliding velocity or maintain surges once started.
Regardless of the model chosen for bed roughness, the necessary conditions for a thick water sheet at temperate glacier beds are:
-
(i) Water-filled subglacial cavities must be rare.
-
(ii) Large quantities of surface melt water must penetrate to the glacier bed.
7. Acknowledgements
I am grateful to Bernard Hallet for guidance and encouragement through the course of this work, as well as for a critical review of the manuscript. Helpful comments on the theoretical analysis were also made by Raymond Fletcher and Paul Delaney.
This work represents part of a Master’s thesis in the Department of Geology, Stanford University.
Financial support was provided by National Science Foundation Grant EAR77-13631 and by the Shell Companies Foundation.
Appendix A. Symbols and Values of Constants
Appendix B. Fluid Flow in the Perturbed Water Sheet
Steady, incompressible flow in the water sheet is described by the equations
where v = 0 has been assumed as has a constant pressure gradient. The velocity w is very small compared with u, but not exactly zero, unless the sheet thickness is constant.
Equations (B–1) and (B–2) can be put in the dimensionless forms
where Ũ and 
are characteristic velocities, l is a characteristic length scale in the x direction, x’ and z’ are dimensionless Cartesian coordinates, u’ and w’ are dimensionless velocities, and other symbols are defined in Appendix A and correspond to usage in the text. Since the dimensionless functions and their derivatives are of O(1). Equation (B–4) demonstrates that 
≈ Ũh/l. Using this estimate of 
, Equation (B-3) can be rewritten as
where (Re) = ρ w Ũh/η w is the Reynolds number for true sheet flow; this is expected to differ very little from (Re) for a slightly perturbed sheet. In Equation (B–5), the left side represents inertial effects, while the term in brackets on the right side represents viscous effects. It is shown in Appendix C that (Re)≤c. 103 for nearly all cases of interest; hence, viscous effects dominate inertial effects if the terms in brackets in Equation (B–5) have a value exceeding c. 104. This condition is satisfied if l/h > c. 104. For h < c. 10 mm, this requires l > c. 0.1 km. a restriction that is virtually always met, since l may be taken as the glacier’s length. Hence, it is justifiable to neglect inertial terms in the equation of motion, which reduces, in dimensional form, to
I now assume that u may be expressed as a mean plus a small perturbation, viz.
where ε⪡ 1. Substituting Equation (B-7) into Equation (B-6) and separating terms of O(1) and O(ε) leads to the two equations:
Equation (B–8a) can be integrated directly:
where à and 
are constants. The perturbation u
1 satisfies the two-dimensional Laplace equation and will have the well-known form,
where the a n , b n , c n , and d n are constants to be evaluated.
The boundary conditions to O(ε) are
The solution for u satisfying these boundary conditions is
Equation (B–12) is used for computing the viscous dissipation Φv (Reference Bird, Bird, Stewart and LightfootBird and others, [c1960], p. 316):
where I have assumed one-dimensional, incompressible flow. From the earlier discussion, ∂u/∂x is assumed to vanish; hence, from Equation (B–12), I can compute
Thus, to O(ε)
Appendix C. Dimensional Analysis of the Thermal Energy Equation
A solution to the complete thermal energy equation can rarely be obtained in closed form. Nontheless, appropriate scaling of the equation can elucidate the important physical processes and show which terms may be neglected.
The complete thermal energy equation for one-dimensional, steady, incompressible flow is (Reference Bird, Bird, Stewart and LightfootBird and others. [c1960], p .315)
where all symbols are defined in Appendix A. This equation can be put into dimensionless form with the scalings
Thus (x’, y’, z’) are dimensionless Cartesian coordinates, u’ is dimensionless velocity, and T’ is dimensionless temperature. This choice of scalings is appropriate as long as “entrance effects”, i.e. very small values of x, are not of interest. Ũ is a characteristic velocity, here chosen to be the mean flow velocity. Substituting Equation (C-2) into Equation (C-1) yields, after rearrangement,
where (Re) = ρ whŨ/η w , (Pr) = η wcw/kw .
The dimensionless derivatives are O(1), so the relative importance of the various terms is determined by the magnitude of the dimensionless constants. Consider first the factor [(Re) (Pr)]−2. (Pr), the Prandtl number, is a material constant equal to 13.7 for water at 0°C, while the Reynolds number (Re) depends upon the flow. In the steady state, melting up-glacier of a point is balanced by discharge, i.e.
Therefore (Re) = ρ w Mx/η w . Table C-I gives (Re) for several values of M and x. The transition to turbulent flow occurs at (Re) = c. 5 x 10 3 (Reference Stuart and Rosenhead.Stuart, [c1963]), so sheet flow ought to be laminar for all cases considered.
Table C-I. Reynolds number (Re) for subglacial sheet flow as a function of melting rate M and distance x along the ice-flow direction
If Re > 0.73, the term [(Re) (Pr)]−2 < 10−2, hence down-stream conduction is quite negligible compared to vertical conduction. Even for 0.23 < (Re) < 0.73, the term is <0.1 and can still be safely neglected. It may be seen from Table C-I that only for very short glaciers and very low melting rates will downstream conduction need to be included in the thermal energy equation.
The factor (kh)2 appears twice in Equation (C-3). Assuming that the wavelength (= 2π/k) of water-sheet perturbations is much greater than h, kh must be much less than unity; hence, terms multiplied by (kh) 2 are negligible in the overall energy balance.
Equation (C-3) can now be interpreted physically. Basically, geothermal heat at the bed and viscous energy losses in the water sheet are available for melting of the basal ice. Part of this energy is not conducted across the film but rather advected down-stream. Corrections for down-stream and lateral conduction are usually extremely small.
Appendix D. Temperature Distribution in the Perturbed Water Sheet
The complete thermal energy equation for one-dimensional, steady, incompressible flow is (Reference Bird, Bird, Stewart and LightfootBird and others, [c 1960], p. 315)
where all symbols are as used in the text and defined in Appendix A, Equation (D-1) can be considerably simplified. Downstream conduction is neglected in accordance with the results of Appendix C. Furthermore, the water temperature is assumed to differ from the pressure-melting temperature by a small enough amount for the following approximations to be valid:
Since P g has been assumed to be constant in the solution of the flow problem (Appendix B), ∂T/∂x is also constant. This is consistent with the neglect of the term in ∂2 T/∂x 2. Finally, the term (∂u/∂y)2 is dropped, since as shown in Appendix B, it is negligible compared to the term (∂u/∂z)2. With these approximation, Equation (D-I) reduces to
The term ∂2 T/∂y 2 has been kept, even though small compared to the other conduction term, because cross-stream conduction might be expected to affect the sheet stability. I assume that the temperature field can be decomposed into a mean and a perturbation, viz.
where ε ⪡ 1. Substituting Equation (D-4) into Equation (D-3), using Appendix B, Equation (B-l1) for the term ∂u/∂z, and separating terms of O(1) and O(ε), the following two equations are found:
The boundary conditions on the temperature are
Since p(x) = p(x = 0) – Pgx, Equation (d-7a) can be rewritten as
where Tmo = –CtP(x = 0) is the pressure-melting temperature at x = 0. Equation (d-5) can be integrated with respect to z to yield
where C 0(x) and C 1(x) are, in general, functions of x to be evaluated by use of boundary conditions.
The temperature perturbation T 1 will have the form
where F(z) and G(z) are chosen to satisfy Equation (D-6). Substituting Equation (D-9) into (D-6), one finds
where f 0 and g 0 are constants to be determined.
The evaluation of C 0(x), C 1(x), f 0, and g 0 is extremely tedious, involving lengthy algebraic manipulations. I present here only their values:
The heat flux at the ice-water interface can now be found by substituting the expressions (D-11) into Equations (D-8) and (D-9), computing the derivatives ∂T 0/∂ z and ∂T 1/∂ Z , and evaluating at z = z0 . The resultant heat flux Qz 0 can be shown to be
If kh ⪡ 1, this can be rewritten by using the expansion of tanh kh for small kh: tanh kh = kh –k 3 h 3 /3 s+ 0(K 5 h 5). With this approximation, the heat flux Q zo can be written in the form of Equation (15) of the main text.
Appendix E. Viscous Relaxation of the Ice-Water Interface
Sagging of the basal ice when the ice-water interface is perturbed is analyzed by the method of Reference FletcherFletcher (1977). Both ice and water are assumed to be incompressible, linear-viscous fluids. The unperturbed state is hydrostatic equilibrium. Imposition of a gentle waviness on the interface results in small changes in the stress state in the fluids, which flow quasi-statically in response to the altered stress condition. The stresses and velocities associated with the perturbation can be written in a form analogous to Fletcher’s equations (17a-d):
where vf and wf are velocities in the y and z directions, respectively, the interface is at z 0 (y) = h + A sin ky, and the slope kA ⪡ 1. A’, B’, C’, D’ are integration constants, and ηf is the fluid viscosity. Equations (E-1) are exact to O(kA). One such set of four equations can be written for each fluid. The boundary conditions used to determine the total of eight integration constants can be stated as follows:
-
(i) Normal and shear stresses and normal velocities in the two fluids must match at z = z0(, Y )
-
(ii) v = W = 0 in the water on z = 0,
-
(iii) v→0, w→ 0 in the ice as z→∞.
The problem of determining the integration constants can be reduced to the solution of a system of four simultaneous algebraic equations. These equations are solved by standard, though tedious techniques.
Due to the extreme viscosity contrast between ice and water, the solution effectively reduces to that for the relaxation of perturbation on a free surface of ice. The quantities of most interest, the velocities, can be expressed to an excellent approximation as
where v p and wp are the velocities in the ice and η i is the effective ice viscosity. On the interface z = z 0(y) these reduce to
where I have expanded the exponentials as power series in kA. Equations (E-3) are the result used in the main text.
