Temperate ice permeability, stability of water veins and percolation of internal meltwater

Abstract

In temperate glacier ice, in situ, besides water veins, there are water lenses, on grain boundaries more or less perpendicular to the direction of maximum pressure p₁ (at the grain scale). Geometry of veins is developed. Grains are modelled as equal tetrakaidecahedra. The stress and temperature fields around a vein at a smaller, microscopic scale are estimated and the water discharge by a Vein is calculated. The time-derivative of the cross-sectional area S of a vein is governed neither by energy dissipation in the water nor by plasticity, but by capillarity effects and salinity. A “vasodilator threshold” p_d for water pressure p_w in the veins is defined. Normally, P_w < P_d, then S has a stable value, the same for any orientation of the vein, and the microscopic temperature is uniform. The coefficient of permeability is proportional to (P_d-p_w)⁻⁴, and thus a true Darcy law does not hold. As an application, the percolation of internal meltwater is studied; in an upper boundary layer about 2 m thick this meltwater flows upwards, because in the bulk of the glacier p_w is very close to P₁, whereas it is zero at the surface. When, exceptionally, p_w > p_d, S increases irreversibly. Whether it leads to the formation of “worm-holes” is discussed.

Numerical Values

• Rheological parameter: B = 0.44 bar⁻³ a⁻¹ = 1.391 × 10⁻²³ Pa⁻³ s⁻¹.

• Lowering of melting point with pressure (air-saturated water): С _m = 9.8 × 10⁻⁸ K Pa⁻¹.

• Lowering of melting point with total curvature of the interface (air-saturated water): C _r = 2.7 × 10⁻⁸ K m.

• Lowering of melting point with salt content: C _s = 1.85 K mole⁻¹ kg.

• Thermal capacity of water: C _w = 4216 J kg ⁻¹ K⁻¹.

• Thermal conductivity of ice: Kⁱ = 2.12 W m ⁻¹ K⁻¹.

• Melting heat: L = 3.35 × 10⁵ J kg⁻¹.

• Interfacial energy between ice and water: γ = 0.034 J m⁻².

• Water viscosity: η_w = 0.0018 Pa s.

• Ice density: ρ_i = 915 kg m⁻³.

1. Introduction

Ice permeability is not a mere curiosity providing an occasion for the advanced student to discover the subtle laws of capillarity. It is a key factor in several glaciological problems.. According to Nye (1976), the outburst of glacier-dammed lakes, in particular the well-known Icelandic jökulhlaups, should proceed from irreversible broadening of the capillary veins. If his theory is correct, this process might explain the disappearance of temporary pools that appear on some glaciers in spring. In any case, the formation of some observed “worm-holes”, several millimetres thick and half a metre long (Raymond and Harrison, 1975), ought to be explained. Also, I have suggested (Lliboutry, 1993) that the permeability of bottom ice replaces the speculated water film between ice and bed, allowing meltwater from the stoss sides of obstacles to reach their lee sides.

Nye’s (1976) theory, which develops Shreve’s (1972) ideas, contradicts Nye and Mae (1972). In Lliboutry (1976) I stressed salinity effects, whereas in Lliboutry (1993) I ignore them. Nye (1991b) examines how salinity modifies the apparent latent heat and diffusivity of temperate ice, neglecting the fact that most of the liquid phase is found in isolated lenses. Many theorists assume that temperate glacier ice has some well-defined coefficient of permeability. Almost none of the authors discuss the dramatic differences between the artificial and stress-free ice they consider and true glacier ice in situ. This confusing literature encourages me to present a clear exposition of the whole topic, instead of limiting my contribution to new ideas and theoretical results.

2. Pressures and Temperatures at the Microscopic Scale

In temperate ice the liquid phase is found as lenses at two-grain boundaries, and as veins at three-grain boundaries (Lliboutry, 1971; Nye and Mae, 1972; Raymond and Harrison, 1975; Mader, 1992a). At any ice water interface, the temperature and the pressure in the liquid phase are well defined. At a curved interface, because of surface energy effects, this water pressure (p _w) differs from the pressure in the ice at the interface (p _a). At a plane interface, the Celsius melting temperature is θ = −C_m P _n, but at a curved interface it is neither −C_m P _w nor −C_m p _n. Let 2/r _L > 0 be the total curvature of spherical lens walls, and −l/r _v < 0 be the total curvature of cylindrical vein walls. The mean salinity of water lenses, and the salinity of connected water veins (in mole equivalents per unit mass, since dissolved salts are dissociated into ions), will be denoted i _L and i _v, respectively. Dissolved carbon dioxide will also be included in the salinity. The physical laws of capillarity read (Lliboutry, 1987, p. 160; Mader, 1992b):

For a lens;

(1)

For a vein:

(2)

The value of the surface energy between ice and water (γ) and of the three constants (C _m, C _r, C _s) entering the formulas are given above. Comparing both Equations (2), P _n, may be eliminated:

(3)

For brevity, a new physical constant will be introduced:

(4)

To go further, it is essential to specify the scale at which Stresses are defined. In the formulas above, p _n is the normal stress at the literally microscopic level (say, at 1 μm from the interface). It may differ from the stress at the grain level, say at 1 mm from the interface, for instance when the vein shrinks plastically. This stress at the grain level, which is not affected by the microscopic processes ongoing close to a vein, may differ in turn from the macroscopic stress, at the decimetre level, because the deformation of a crystal is very anisotropic, and identical strain rates and stresses in all the individual crystals are impossible.

When studying strain rates under an applied stress, or the reverse, on polycrystalline samples (such as ice cores), it is the stresses and strains at the 0.1 m scale that are considered. For modelling glacier flow, a creep law at a still larger scale, say 10 m, is needed, because in a glacier there are “blue bands”, fine-grained layers, and overall fluctuations in the water content. I shall as usual not distinguish between stresses at the 10 m and at the 0.1 m scale, or between the latter and stresses at the grain scale. The reasons for assuming that the stress is the same in every crystal (whereas strain rates differ) are given in (Lliboutry 1987, p. 456; 1993, p. 53). The stress down to the millimetric scale will be called the local stress, whereas the stress near veins or lenses, which is of interest here, will be called the microscopic stress.

In unstressed ice, p _n is roughly the same everywhere and in any direction. When i _L = i _v, both the water pressure and the temperature are lower in the veins than in the lenses. Thus, heat goes from the lenses to the veins. Countinuously water freezes at the walls of lenses and the walls of veins melt, until the former disappear. The microscopic stresses due to the Volume increase when a lens freezes are fast relaxed (Nye and Mae, 1972).

(The reader who finds it paradoxical that the melting temperature is lower in the veins, where water pressure is lower, should consider the case of a lens connected with a vein. If the temperature were higher in the vein, heat would diffuse steadily through the ice from the vein to the lens. There would be melting at the lens walls, meltwater flow into the vein, and refreezing of the same mass of water at the vein walls; and so on indefinitely, without the need to supply energy to the system. Such a perpetual process would be contrary to the second principle of thermodynamics.)

When a sample of polycrystalline ice is kept for a long time in a water bath (or left in a cold place slightly above melting temperature), the surface of the sample acts as heat source, and the veins go on coarsening (Nye, 1991a: Mader. 1992b).

Within a glacier, the slate of stress is no longer hydrostatic, and thus p _n may vary from one water-ice interface to another. For water lenses at a two-grain boundary, p _u is highest and θ _L Lowest when the boundary is more or less perpendicular to the direction of maximum normal pressure. These lenses enlarge at the expense of others that freeze, and thus only lenses with this orientation will be considered. Denoting the three principal local stresses, with reversed sign, p ₁ > p ₂ > p ₃, then, shortly after the final size has been reached, p ₁ = p _n, and the walls of these lenses are at temperature:

(5)

(There is a printing error of sign in Lliboutry (1993, equation (6)).) The “locally stress-controlled temperature” that I consider in this paper does not differentiate between θ _v and θ _L, because the influence of salinity and capillary effects could he neglected for that study. Now, to study permeability, more refinements, and the temperature field at a smaller, truly microscopic scale, must be considered.

(By the way, note that in a temperate glacier the local ice temperature decreases with depth by about C _m ρ _i g m⁻¹. Thus a very small heat flux goes down, which causes an imperceptible melting rate K _i C _mρ_i g/ρ_iL) = 0.192 mm a⁻¹ at the sole.)

A vein has three cylindrical sides (its “walls”), each one limiting a different crystal. (A vein with four sides has a vanishing probability of being found. For this reason an assemblage of equal cubes would be a very bad model for polycrystalline ice.) In a steady regime, the temperature at the three walls must be the same, ρ _v. According to Equation (3), the curvature of the three walls must be the same, 1/r _v. Therefore, from Equation (2), the microscopic pressure within the ice at the walls is the same. p _n, whereas at a large distance the stress is anisotropic.

Locally, since the veins are connected together, p _w and i _v are the same for all veins (or everywhere), and thus ρ _v is lowest where r _v is smallest. This means that veins with the smallest r _v catch more heat from the surrounding ice than do others. Therefore, if all the cross-sections of veins were similar (as will he assumed later), all the veins Should tend to have the same size locally, whatever their orientation, in contrast to lenses.

3. Geometrical Relations

In fact this similarity is not observed, because the surface energy between two crystals varies with the relative crystal orientation (Walford and Nye, 1991). It is lower when coincidence site lattices form, and the multimaxima fabric that is observed in temperate glacier ice should come from this formation (Higashi, 1978; Matsuda and Wakahama, 1978).

Consider three crystals, 1, 2, 3, having an edge in common, with “grain dihedral angles” at this edge w ₁, w ₂, w ₃, respectively (Fig. 1a). The surface energy of the boundary between crystals i and j will be denoted γ_ij. If the grain dihedral angles are governed by surface energies only, the three vectors in the directions of the three grain boundaries, and with respective intensities γ₁₂, γ₂₃, γ₃₁, must have а zеrо vectorial sum (Fig. 1b). Consequently:

We must have the inequalities:

(6)

and two similar ones obtained by cyclic permutation of the indices.

The sines relation in a triangle yields:

(7)

From the cosine relation in a triangle:

(8)

and two similar ones.

Fig. 1. Grain dihedral angles at a three-grain intersection, as governed by interfacial energies.

These relations are valid in stress-free annealed ice. In a glacier, the self-energy of dislocations perturbs the equilibrium of the interfacial forces, grain boundaries migrate, and ice recrystallizes continuously. For the time being, due to a lack of studies on quenched samples of glacier ice, it will be assumed that Relations (7)–(8) are also valid in situ.

When a vein with cross-section ABC exists at the edge considered (Fig. 2), the equilibrium of the interfacial forces must be obeyed at the three vertices A, B, C. Consequently:

Since the interfacial energy between ice and water is independent of the orientation of the crystalline lattice at the interface, the grain boundaries bisect the dihedral angles at A, B, C.

With

denoting the values of these “vein dihedral angles”:

(9)

Fig. 2. Vein cross-section ABC, and centres of curvature of its three walls C ₁, C ₂, C ₃.

The centres of curvature of the three walls are denoted C _l, C ₂, C ₃, and the three-grain intersection without vein is O. Since γ_v = CC₁ = CC₂, the triangles OCC₁ and OCC₂ are equal, and thus ОС₁ = 0C₂ = R. Similarly, OC₃ = OC₂ = R. The angles of triangle C₁C₂C₃ are equal to π-w ₁, π-w ₂, and π-w ₃, respectively. The sines relation reads:

(10)

Consequently, the angle of ОС₁ with boundary (1–2) is π ·w ₃, and

. This is the maximum value that

can have (then the distance OC vanishes). Thus we must have:

(11)

Similar conditions holds for w ₁ and w ₂. Comparing with Equation (7), a necessary condition for the existence of a vein is:

(12)

This condition is not always met. Moreover, it is not sufficient. Another condition is that O be at the interior of triangle C₁C₂C₃. It requires w ₁, w ₂, w ₃ > π/2, or, geometrically, that the triangle with sides γ_ij (Fig. 1b) have three acute angles. Figure 3 shows a low-energy boundary (1–2), and four others with higher energies: about the same for (2–3) and (3–1), different for (2–4) and (4-1). A vein exists at the intersection (1–2–3), but does not exist at (1–2–4) because a grain dihedral angle is acute.

The existence of about 5% of three-grain intersections without veins was discovered by Raymond and Harrison (1975) and confirmed by Mader (1992a). My Figure (3) may be compared With the latter’s Figure (1): a quadrilateral grain boundary with a measured low energy had veins on two sides, and none on the other two.

Fig. 3. A vein exists at the triple intersection (1–2–3), but not at (1–2–4), because grain 1 has an acute grain dihedral angle.

All the theory above is consistent with Mader’s observations, although the conditions for the existence of a vein differ from hers. According to Mader, it is

. This condition expresses only the fact that the three walls are concave (when seen from outside). They must be so in stress-free ice with a low salinity, for the veins to be at a lower temperature than the lenses. But if

, the walls might be convex. It will be shown that no stable value of the cross-sectional area exists in this case, however.

To develop the theory taking into account different interfacial energies between neighbouring crystals would be extremely difficult. Therefore, the study of the behaviour of a vein will be done assuming that the three γ_ij are equal, with the hope that in general they are not very different. The cross-section of the vein is then an equilateral curvilinear triangle.

The problem is to decide the best numerical value for this grain-boundary energy, or, according to Equation (9), the best value for the dihedral angle Ψ. The first published values span 20–34°. An improved analysis of Walford’s data by Nye yielded Ψ = 25° (Walford and Nye, 1991). On the other hand, Mader (1992a) found once, in the quoted case of two veins out of four missing at the periphery, Ψ = 105°. (She used then spikes entering at its extremities a three-grain intersection without vein.) All these measurements were done on annealed ice, produced artificially by a technique that should favour the formation of a fabric but not the one found in temperate glaciers. The value Ψ = 30° has been adopted in the following numerical calculations.

The cross-sectional area is readily found to be (Nye, 1976):

(13)

It is also of interest to know the ratio of the radius r _i that inscribes the cross-section to the radius r _c that circumscribes it. It is found:

(14)

Some values are given below:

4. Water Discharge by a Vein

The abscissa along a vein is denoted s, and measured iu the upstream direction. Let Z be the altitude of a cross-section of the vein (relative to some reference level,). The hydraulic head of the vein at that point is defined (ρ_w denoting water density and g gravity) as:

(15)

The slope of the vein (positive when water flows downwards) is denoted tan β. The hydraulic gradient reads:

(16)

Flow in a vein is obviously laminar. Then, from Newtonian viscosity theory, in a transversal plane, η_w denoting water viscosity, the velocity of water (u) obeys:

(17)

At the periphery of the conduit, u = 0. Thus. u is proportional to A, and so is the discharge Q. For dimensional reasons,

(18)

where u denotes some numerical coefficient.

For a circular cross-section whose equation in polar coordinates is r = a, we have the well-known Poiseuille solution, which is used by Nye (1976):

(19)

There is another simple solution of Equation (17) when the cross-section is an equilateral triangle, with sides equal to 2α and area

. Adopting Cartesian coordinates such that the three summits are (−a, 0), (a, 0) and (

), the polynomial of lowest degree that is zеrо along the triangle is:

(20)

For this special triangle the Laplacian reduces to a Constant,

, and thus C may be chosen so that Equation (17) is satisfied at any point. It is found:

(21)

The ratio r _i/r _e is 1 for a circle, 2 for an equilateral triangle. Both are special cases of curvilinear equilateral triangles. For all of them, intuitively, μ is a continuous function of r _i/r _e. Therefore, although a precise numerical Computation of μ for a vein might be done, the value obtained by a mere linear extrapolation, namely μ = 0.0264, has been confidently adopted.

5. Evolution Equation for a Vein

The variation of S with time may be written:

(22)

where E is the heat produced by viscous dissipation in the water,

is the heat received by the walls of the vein from the surrounding ice (both per unit length and unit time), ρ_i L is the melting heat per unit volume, and P is the rate of shrinking by plastic deformation.

The first term has been calculated by Röthlisberger (1972). He took into account the fact that part of this heat is used to warm the water, when it flows towards lower pressures and higher melting temperatures. His value for the dimensionless constant C _m C _wρ_w that gives the ratio has been modified (Lliboutry, 1983), because C _m is larger when water is air-saturated, as should be the case given the air inclusions in glacier ice (Raymond and Harrison, 1975). This ratio is 0.413 instead of 0.316. Nye (1976) noted that the water-warming rate involved is not

, but the material derivative

. The corrective term is totally negligible, however. Taking into account Equations (16) and (18), we obtain:

(23)

The term P comes from the fact that at a vein wall the microscopic normal stress is p _u, whereas at some distance it is the local stress. The fact that each of the three crystals surrounding a vein deforms differently, according to the orientation of its lattice, makes illusory the classical formula used for Röthlisberger channels. I shall use it, however, to obtain an order of magnitude. It will show that this term is negligible.

For this reason, and because surface energies, which act to preserve the shape of the cross-section, were not considered, a smart calculation by Oakberg (1981) is not very meaningful. (He found radial velocities at the middle of the sides of the cross-section, v(r _i), that are 2.46 times the velocities at the corners, v(r _c). For the cross-section to keep a constant shape, the ratio must be r _i/r _c = 0.397 instead.)

The component of the local stress that intervenes is the radial pressure p _r in the plane of the cross-section, say the xy-plane. It varies with azimuth ϕ as:

(24)

As microscopic pressure at a large distance, the mean value Of

, is adopted. Let α, β, γ be the cosines of the Z axis (i.e. of the vein) with respect to the principal directions of stress. Then:

(25)

The microscopic pressure at the walls of the vein is p _n, as given by Equation (2). The classical formula for the closing of a cylindrical conduit, when the medium has a third-power-law viscosity, namely:

(26)

reads:

(27)

Lastly, the heat flux reaching the vein must be estimated. Grains are modelled as equal tetrakaidecahedra, each one formed by eight regular hexagous and six squares, all with equal sides (the 36 edges of the polyhedron). This model was used by Kelvin in 1867 in his speculations about the structure of ether, and in modern times by Nye and Frank (1973). It is displayed in Figures (4) and (5). If b is the length of the edges, any edge is at distance

from the centres of its two neighbouring hexagons, and at distance b/2 from the centre of its neighbouring square. The arithmetic mean of grain areas, as measured on sections, is

when the sections are parallel to hexagons, and

when they are parallel to squares. Thus, the distance from a vein to the centre of a neighbouring grain boundary is

and

, respectively.

Fig. 4. Tetrakaidecahedron obtained by truncating the six vertices of a regular octahedron with edges equal to 3b, so that the central thirds are left. It is formed by eight regular hexagons and six squares, and has 36 edges of equal length b. Thе octahedron has a volume

, and each truncation removes 1/54 of its volume. Thus the volume of the tetrakaidecahedron is

. Its surface has a total area

, and S³/V² = 150.12. (This ratio is 216 for a cube, and 113.10 for a Sphere.) (a) Plan view, and (b) elevation when a square is horizontal. The height of the polyhedron is then

, and the mean area of horizontal sections is 4b². Edges that are not horizontal are 45° from horizontal. (c) Plan view, and (d) elevation when a hexagon is horizontal. The height of the polyhedron is then

, and the mean area of horizontal sections is

.

Fig. 5. Infilling of space by a regular assemblage of equal letrakaidecahedra. A single layer is represented. The layers below and above are in contact at the square left white. The grain dihedral angle between two hexagons is

; that between a square and a hexagon is

. For Equation (8) to be obeyed, the ratio of the surface energies between two square boundaries and between two hexagonal boundaries must be

.

I assume that the heat flow is the same as if the temperature were uniformly θ _L at disiance

, and θ _v at disiance

from the vein axis. The sought heat flux is then

(28)

where K _j denotes ice thermal conductivity, and θ _L and θ _v are given by Equations (5) and (3), respectively. Since

, the term 2C_R/r _L will be neglected. As for r _v, it is given by Equation (13).

Equation (22) can now be written in full, as:

(29)

Numerically, with the values given at the beginning of this paper, μ = 0.0264, and V = 0.284:

(30)

To suppress many powers of ten, S may be expressed in units S₀ = 10⁻¹⁰ m², p ₁, p _v, and p _w in bar = 10⁵ Pa, and i _L, i _s in μmole g⁻¹ = 10⁻³ mole kg_-1. These units are in the same order of magnitude as the corresponding variables. With them, Equation (29) reads:

(31)

6. Discussion

With the values of the hydraulic gradient and of p _v - P _w that may he expected on any glacier, the terms E and P of Equation (22) are normally totally negligible, contrary to Shreve’s (1972) opinion, and to Nye’s (1976) theory. With only these terms, the steady state

was unstable against any perturbation in S. With the predominan terms coming from

(i.e. from capillary effects), in contrast, the steady state is stable: if S increases, S ^−1/2 decreases and

becomes negative. (Note that the existence of a Steady value comes from the fact that the vein walls are concave. If they were convex, Equation (13) would remain valid, but in Equation (29) v would be replaced by −v. Then

would be negative for any S.)

Define:

(32)

The steady state corresponds to a uniform temperature at the microscopic scale. The steady value ofS is:

(33)

It is independent of the orientation of the vein, as already stated. Nevertheless, this steady value of S is stable against a perturbation in the water pressure only if P _w < p _d. Then, if P _w increases, S increases until

vanishes again, and the steady regime is restored. This Stability ends, and the steady value of S becomes infinite, when p _w becomes equal to p _d. Since blood vessel was the original meaning of “vein” and humour is not forbidden, p _d will henceforth be called the “vasodilator threshold”.

It is interesting to examine whether a very strong hydraulic gradient, such as is artificially created when a water-filled borehole reaches the vicinity of a subglacial waterway or cavity at atmospheric pressure, modifies p _d significantly. At Glacier d’Argentière, French Alps, a borehole emptied when reaching 120 m in depth, and it remained empty when the bed was reached, 5 m below (Hantz and Lliboutry, 1983). Thus the hydraulic gradient was about 20. Equation E in Equation (22), but always neglecting P, Equation (1) reads:

(34)

where

is the value of p _d above. Since

, sin β may be neglected, and:

(35)

With the term E included, the function f(S) has now a minimum for S = S _m. In general, f(S) = O has two roots, and the smaller one is the stable solution above. When f(S _m) > 0, there are no roots, and S grows continuously, at an increasing rate. Since σ/S is very large, λ is a very mild function of S, which may be considered as a constant when calculating the minimum S _m. (Next, for obtaining numerical values, the precise λ corresponding to each S _in, say λ_m, as given by Equation (28) will be used.) The result is:

(36)

The vasodilator threshold is the value of p _w that makes f(S _m) = 0:

(37)

Numerical values of the corrective term p _d - p _d°, for

are given below:

The corrective term is negligible for any hydraulic gradient.

Plausible values of (i _L — i _v) will now be assessed.

At the atomic scale, internal meltwater appears when a dislocation and its elastic self-energy disappear, a process that occurs mainly at grain boundaries. Since the impurity concentration is larger at the grain boundaries, internal meltwater should be less pure than surface meltwater coming from the total melting of snow and ice.

First, consider measurements on unstressed ice cores, such as those by Raymond and Harrison (1975). In this case (p ₁ — p _w) was zero. Since thermal boring and subsequent melting at the surface of the core should have flushed out the salts from the veins, i _v = 0. Then, from Equation (32), i _L should be related to the observed cross-section by:

(38)

S was found to be about 7S ₀, and thus i _L = 0.176 μmole g⁻¹, If w is the water content of the ice core, almost all due to the water lenses, the global salinity of the ice core is

. From the electrical conductivity of the melted cores, Harrison anil Raymond (1976) estimated the global salinity as 3.3 × 10⁻³ μmole g⁻¹ for fine ice (grain diameter= 2 mm), 0.87 × 10⁻³ μmole g⁻¹ for coarse ice (10–50mm). These values yield w = 1.88% and 0.49%, respectively. They are quite acceptable, but, since the water content was not measured, the present theory cannot be tested.

The study above was done on ice cores from Blue Glacier, in the Olympic National Park, Washington, U.S.A., a very wet area. In the Alps, where precipitation is lower, the salt content of snow and glacier ice is higher by a full order of magnitude (Lliboutry, 1976).

Nevertheless, in situ. i _v is no longer negligible. Veins are fed only in part by surface meltwater, especially since shallow glacier ice is very bubbly, and air bubbles clog water flow in the veins. A large part of the water feeding the veins is internal meltwater. In stressed ice, water lenses at grain boundaries must grow continuously until they reach a vein and can empty into it (Lliboutry, 1987, p. 123). In the extreme case of veins exclusively fed by internal meltwater (a model adopted by Berner and others (1978)), i _v = i _L. Then, the vasodilator threshold depends on the maximum compressive stress p ₁ only.

When p _w exceeds the vasodilator threshold, S grows continuously. For S in the order of 1000

, S^−1/2 has practically vanished, and, even for moderate values of the hydraulic gradient, in Equation (22) the term E must be taken into account. It affords a positive feedback. Consider several distinct paths between two given four-grain intersections. The variation of the hydraulic head Z_w is the same for all paths, and thus

is lower the more tortuous the path. Consequently, the straightest path thickens fastest, and the feedback along it is the most important. The global water flow is no longer equally shared between veins of the same cross-section. It tends to concentrate into several almost straight water-channels of millimetric thickness. Such ducts were observed by Raymond and Harrison (1975). I call them “worm-holes”.

It remains to assess how long it takes for such worm-holes to form, by integrating:

(39)

assuming the realistic values

, sin β = 1, and E ₁ = 1.8 × 10¹¹.

With the first term alone, a worm-hole should take several centuries to appear. With the term in (p _w—P _d) alone, S might reach 1 mm² in 1 d, but S would grow equally for every vein, and worm-holes would not form. Therefore. I suggest a two-step model. First, р _w > P _d, and all the veins thicken. Next, p _w = p _d, and the worm-holes form. There are not enough observations to clear up definitively this puzzling problem.

In the quoted case of an extremely large hydraulic gradient at the bottom of a water-filled borehole, the formation of worm-holes through which it emptied should be expected. It seems a more credible process than any other hypothesis that comes to mind: thickening of the veins by the influx of warm water, hydraulic fracture, or fortuitous formation of a fault there.

6. Ice Permeability

In polycrystalline temperate ice all the veins are Interconnected, when they are not clogged by air bubbles. Thus, at the local scale, an interstitial water pressure p _w and the corresponding hydraulic head Z _w and hydraulic gradient

can be defined. As long as the vasodilator threshold is not reached, the total water flow per unit area (q) may be expressed as in Darcy’s law:

(40)

but the coefficient of permeability k does not depend on geometrical factors only, as in soils or porous rocks. It depends also on the local values of p ₁, p _w and (i _l-i _v).

To pass from discharge by a single vein Q to the water flux |q|, some geometrical model for the polycrysial is needed. It will be assumed that the grains are tetrakaidecahedra of equal size. It may be objected that they are not the polyhedra that ensure a minimum surface for a given crystal volume (Weaire and Phelan, 1994). This objection has little weight. It is not the area of the surface that is minimized by capillary forces, it is the total surface energy. Equation (8) is obeyed if the surface energies of square and hexagonal boundaries, say

and

, are in the ratio

. Note ihat the orientation of the hexagons, which are the low-energy boundaries, recalls the four-maxima fabric.

The objection that grains in polycrystalline ice are of different sizes might give more concern. With equal tetrakaidecahedra, water contents of 1% or more would be impossible, because the water lenses would then be larger than the square sides. I believe that this geometrical model leads to good estimations, however, when statistical parameters other than the water content are sought.

The discharge by a single vein, when S has its steady

(41)

To obtatn the analogue of Darcy’s law, this value must be multiplied by three factors:

• 1. A vein density factor. The number of veins per unit area that cross a section of polycrystalline ice increases with the number of grains counted in this section. With σ denoting the mean cross-sectional area of grains, this number of veins (the “vein density”) is proportional to σ^−l. With the model of equal tetrakaidecahedra, a vein density of 2/σ is found. Raymond and Harrison (1975), who examined six thin sections, found vein densities ranging from 0.8/σ to 2.1/σ, but they used an arithmetic mean for σ. A harmonic mean would have been more appropriate, and it is smaller.

• 2. A tortuosity factor, because there is never a set of successive veins exactly in the direction of the hydraulic gradient. With the tetrakaidecahedra oriented as in Figure (4)a and b, and in Figure (5). when the gradient of Z_w is vertical,

  
  for all the veins. This value will be adopted. A cumbersome statistical calculation for all the orientations would not be very meaningful.

• 3. A path continuity factor. A series of grain edges is not necessarily a possible path for water flow, for two reasons: at some edges no vein exists, and even if a vein exists, it may be clogged by an air bubble (Lliboutry, 1971). We have not enough data to assess the probability of the first circumstance in temperate glacier ice. We assume that it is negligible. The probability of the second circumstance is easier to handle.

To obtain, among the number of veins crossing a plane, the proportion of those that permit water flow, Raymond and Harrison (1975) considered the case of several bubbles on the same edge. Their result does not differ significantly from that obtained when this case is ignored. If b is the mean length of an edge, and s the mean distance of bubbies along a path (their notation), this ratio is then (1—b/s).

We may go further, by introducing the number of air bubbles per unit volume (N), and their mean diameter (δ). With ρ denoting the density of bubbly ice (860–900 kg m⁻³), and ρ the density of bubble-free ice (915 kg m ⁻³), we have:

(42)

A bubble intersects a path when its centre is less than

from it. Thus the number of bubbles per unit length of path is:

(43)

As for the tetrakaidecahedron model yields

. Consequently, the probability of free passage along a vein is:

(44)

For a given air content, the probability decreases for small bubbles, and increases for fine ice (not as dramatically as the vein density factor, however).

With Q given by Equation (10) (where μ = 0.067), S given by Equation (24) and K given by Equation (30), once the three factors above have been taken into account, the permeability reads, reverting to SI units:

(45)

Mader (1992a) notes that, when the probability is low, there are networks of uninterrupted veins that do not extend to infinity and thus cannot ensure macroscopic permeability. Compulations by Frisch and others (1962) and by Sykes and Essam (1964) have shown that this begins to happen when the probability is lower than 0.6. When it is lower than 0.39 the material becomes macroscopically impermeable.

Thus, Equation (45) should be valid for (l—b/s) > 0.6 only. For the realistic value ρ = 0.87, the condition reads

, which is generally the case.

7. Percolation of Internal Meltwater

The theory above will be used to handle the percolation of internal meltwater. The hydraulic gradient

is assumed to be vertical. Wiih the z-axis vertical upwards, its intensity is then:

(46)

The rate of internal melting equals the increase of volume discharge Q(Z) downwards. They depend on the effective shear stress т. as defined in Equation (26):

(47)

The origin of the zaxis is at the sole, and the surface is at z = h. It is assumed that: (1) no external meltwater enters the vein system at the surface, but instead q(h) is negative: thus i_v = i_L. and p _d = p _w the vasodilator threshold is never reached: thus p = p ₁ — p _w is every-where positive. Both assumptions will be checked once the solution has been found.

It is also assumed that p ₁ = ρg(h — z) + T, as in the plane problem, and that T is a constant for any z. In fact T is in the rangе 0.4–0.7 at the surface, and equals about 1 bar near the bottom, but to assume a constant value does not change qualitatively the results, and makes the calculation simpler and its difficulties clearer.

With Z _D denoting the unknown level where q = 0 (D is for “water Divide”), Equation (47) is readily integrated:

(48)

and p reads:

(49)

On the other hand. Equation (40) reads:

(50)

From Equation (46) and (49), p _w and Γ are related to dp/dz by:

(51)

Comparing with Equations (48) and (50):

(52)

Function p(z), as q(z) and P _w (z), obeys a differential equation of second order, subject to two boundary conditions: P _w (h) = 0 and p _w(0) = pgh. Equation (52) is a first integral, including the integration constant z _D. It is convenient to adopt a unit pressure

and a unit length Δ such that:

(53)

Adopting ρ/ρ_w = 0.9, k ₁ = 0.02, and T = 10⁵ Pa, it is found:

(54)

The relation between the dimensionless variables

and Z = z/Δ reads:

(55)

P has been assumed to be positive for any Z, and must remain finite, as is P ₁. At the boundaries Z = O and Z = h/Δ = H, P is subject to conditions

(56)

A value H = 60 (h = 294 m) has been adopted tor the numerical computation.

The fact that the differential equation is strongly non-linear makes it impossible to solve by the usual methods. I have found the following one, but I leave more competent mathematicians to assess its accuracy.

The non-linearity leads to the formation of boundary layers where |dP/dZ| is very large, whereas in the bulk of the glacier P and dP/dZ are smaller than 1. In all this central part, P can be calculated as a function of Z — Z _D = X, without the need of some given value of P to start with, thanks to the following algorithm:

(57)

This algorithm converges very fast towards a solution P(X) that does not depend on the the values of P in the boundary layers. P increases slowly from 0.36541 for X = −56 to 0.69966 for X = −4, and dP/dX increases slowly from 0.00163 to 0.04147 in the same interval. Next, starting from this solution, P(X) in the upper boundary layer can be determined

(58)

This equation has been solved by a Runge-Kutta method. It is found that P = 20.8 is reached when X = 0.4001. Since then dP/dX = 7,3860, this value of X is quite insensitive to the precise value of P at the surface. Thus. H-Z_D = 0.4001 (= 1.96 m). At Z = Z _D, P = 1..3562 (= 0.065 bar), and dP/dZ = 1. 1 shall define a boundary layer as a region where |dP/dZ| > 1. With this definition, the upper boundary layer is the one whose internal meltwater exudes at the surface. (It exudes 0.26 mm a ⁻¹ with the adopted value T = 1 bar, but only 0.016 mm a ⁻¹ if T = 0.5 bar.) No external meltwater enters the vein system and lowers its salinity.

Now. as Z _D has been determined. P(Z) can he calculated in the lower boundary layer. Equation (55) reads:

(59)

From Z = O to Z = 0.004, since dP/dZ is very negative, the following approximation is accurate enough:

(60)

Here also, the precise value of P at the boundary does not significantly modify P(Z).

Next, starting from Z = 0.004. P = 1.11815, the Runge-Kutta method has been used. It is found that dP/dZ = −1 at Z = 0.094. The thickness of the bottom boundary layer (as I define it) is thus 0.094 = 0.46m. Next, at Z = 0.65. dP/dZ is zero, and P has its minimum value, 0.3609. (The corresponding maximum value of the vein’s cross-sectional area is S = 0.06615 mm².) The values of P given by Algorithm (57) are found again around Z = 2. Thus, for Z < 2, we have an interesting example of an algorithm that converges very fast towards a wrong solution.

8. Conclusion

It has been shown that capillary effects and salinity, which exclusively control the size of veins (contrary to the case of Röthlisberger channels), in general ensure a steady Value. Nevertheless, since this value depends strongly on the water head in the veins, hydraulics are quite different from those governed by Darcy’s law. Given the discharge (for example, when percolation of internal meltwater is considered), it is much more the water head Z _w and the corresponding permeability, rather than the hydraulic gradient, that adjust to ensure the discharge. With Darcy’s law instead, it would be the hydraulic gradient only. The water pressure p _w differs from the maximum compressive stress by no more than 0.065 bar. Consequently, for the boundary conditions on p _w to be obeyed, boundary layers with large hydraulic gradients, but very small vein sizes and permeabilities, must appear.

This unexpected behaviour, of a material that becomes more impermeable at its surface to maintain a high internal pressure, would make the fortune of a tyre manufacturer. For a glaciologist trying to determine the permeability of temperate ice, it is a disaster. On a stress-free sample, permeability is totally different from in situ. In an in situ experiment, since a very high hydraulic gradient must be applied in order to assess the permeability within a reasonable time, the vasodilator thresh-old may be exceeded, and the result may be erroneous. Probably for this reason, old experiments in situ of the 19th century (quoted in Lliboutry, 1971) lead to the idea that clear temperate glacier ice is quite permeable.

In fact, the in situ permeability of temperate glacier ice is very small indeed, except when it is bubble-free and with very small grains, or when internal meltwater production is very large. Such conditions are met at the sole of glaciers sliding fast on a rough hard bed. Therefore, my theory about internal melting and migration of the meltwater that allows a refreezing trend at the bed (Lliboutry, 1993) is not contradicted by the present theory.

The layers of clear ice that form the popular “blue bands” at the glacier surface may also have an appreciable permeability. These layers come from the concentration of the strain by a feedback process that provides small grains. Nevertheless, when they are a relict form, recrystallization has enlarged the grain-size, the multimaxima fabric is found again, and only a flattening of the crystals gives evidence of their origin (Vallon, 1967: Lliboutry, 1987, p. 122–123). So we cannot assert that blue bands are paths for water percolation, unless observations in tunnels discover this phenomenon.

Putting aside (1) the obvious case of surface ice after days of fine weather, when individual grains have loosened from each other because of solar radiation, (2) the poorly documented cases of “ice mylonites” along internal faults, and (3) the mysterious “ice worm-holes”, there is no more water percolation through the bulk of a temperate glacier than through silty clay or granite (Brace, 1984). In contrast to Shreve’s (1972) views, surface meltwater reaches the bed only by moulins, old closed marginal crevasses and other large waterways.

Another consequence of the present theory is that water pressure in subglacial cavities has nothing to do with the interstitial pressure in the vein system.

An obvious improvement of the theory would be to use empirical statistical values instead of Values drawn from a simplified model. We need statistics not only of in situ values of the vein dihedral angles Ψ, but also of the corresponding grain dihedral angles w. This is a challenge to the ingenuity of experimental investigators. Nevertheless, the best test of the theory would be to measure in situ the interstitial water pressure p _w without disturbing it. This may be done by putting pressure gauges at the bottom of boreholes, and carefully sealing them with cores of cold ice.