Skip to main content Accessibility help


  • Access


      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        The Flow of Ice, Treated as a Newtonian Viscous Liquid, around a Cylindrical Obstacle Near the Bed of a Glacier
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        The Flow of Ice, Treated as a Newtonian Viscous Liquid, around a Cylindrical Obstacle Near the Bed of a Glacier
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        The Flow of Ice, Treated as a Newtonian Viscous Liquid, around a Cylindrical Obstacle Near the Bed of a Glacier
        Available formats
Export citation


This paper describes an analytical solution of the equations of motion and heat conduction for ice flowing around a cylindrical solid inclusion and over a solid plane boundary. This is intended to be a simplified representation of the flow of clean glacier ice around a stone and over a rigid rock bed. The ice is treated as a Newtonian viscous liquid and the equations are solved in two dimensions. Regelation boundary conditions are applied at both ice–rock interfaces. It is found that finite solutions for the temperature and stream function only exist for the special cases in which two dimensionless critical wavelengths are zero. That is, unless the stone is very far from the glacial bed, the classical regelation boundary conditions cannot be obeyed over the whole of its surface.


The problem of describing the motion of a clast in the basal ice of a glacier lies at the heart of physically-based models of erosion and deposition. Order-of-magnitude calculations have been made by Lewis (1960), who discussed the maximum force that could be exerted by a rock held in ice, and by Boulton ([c1975]), who adapted Weertman’s (1957) sliding theory to produce a condition for the onset of abrasion. Röthlisberger (1968) identified and described the processes which tend to bring clasts into contact with the bed. Full analytical solutions for the motion of isolated clasts in ice unbounded by other solid surfaces have been obtained by Glen and others (1957), to explain the preferred orientation of the long axes of clasts in the direction of flow of a glacier, and by Watts (unpublished), in order to derive equivalent viscosities for dilute suspensions of particles in ice. This paper describes an analytical solution for the flow of ice past a solid object and over an external solid boundary. The analysis indicates that interaction between processes occurring on the glacier bed and on the surface of a clast in the basal ice will produce effects which have not been considered in previous calculations.

Figure 1 shows the simple two-dimensional problem that will be analysed in this paper. A long cylindrical stone with circular cross-section, radius R, lies in the ice at the point x = 0, y = ∆. The stone lies in an area zone (i) where the glacier bed and the surface of the ice are horizontal planes, y = 0 and y = H. Up- and down-stream, in zones (ii) and (iii), the glacier bed is undulating and has an overall slope. We shall suppose that the boundaries of zone (i) are far enough from the origin for the ice flow not to be affected by the presence of the stone.

Fig. 1. A cylindrical stone radius R lying a distance ∆ from a plane rock surface in ice of thickness H.

In basal-sliding theories dealing with clean ice (Weertman, 1957; Lliboutry, 1968; Nye, 1969, 1970; Kamb, 1970) it is usual to assume that a lubricating “regelation layer” of water exists between ice and bed-rock. Thus, locally, the shear stress at the ice–rock interface is zero. However, if the bed undulates about a given base plane the ice exerts a net force on the rock which can be interpreted as the result of an average shear stress acting over the base plane. In zones (ii) and (iii) there is an average shear stress between the ice and the undulating bed. In zone (i) we will suppose that the bed can be given an effective or “magic” roughness for part or all of its area, despite the fact that geometrically it is treated as a flat plane. This means that the shear stress on this part of the boundary y = 0 can take non-zero values. Thus, the transition between the region of flow over an effectively rough bed and the region of flow around the smooth stone and over a smooth bed is independent of the choice of the extent of zone (i). This allows a physically reasonable description of the boundary conditions on y = 0 to be combined with simple boundary conditions describing the velocity distribution in the ice at the boundaries between zones (i), (ii), and (iii). The geometry of this problem has been specified in such a way that the two boundaries at which the physical processes are most complicated, the ice–rock interfaces, are coordinate curves in the bipolar coordinate system.

Suppose that in zones (ii) and (iii) the ice moves under gravity with a steady-state flow pattern which depends on the shape of the bed and the depth of the ice. In general there will not be a steady-state situation in zone (i), for, unless the stone is at rest with respect to the bed or moves in such a way that the separation ∆ is maintained, the pattern of flow must change with the change in the relative positions of the ice–rock interfaces. This problem does not arise in the analysis of the motion of isolated clasts in ice.

For the sake of making a relatively simple first analysis, it is assumed in this paper that forces — X, —Y and a couple L oppose the forces X, Y, and couple L exerted by the ice so that a steady-state situation is maintained. Thus, only a small sub-set of the possible solutions for the motion of the ice and stone is being considered. The set of solutions with X = 0, L = 0, and — Y equal to the weight of the stone is physically possible. Other solutions are artificial, but the direction of the resultant force that would have to be applied to the stone, in order to maintain its position with respect to the bed, gives an indication of the direction in which the stone would tend to move in the real situation where gravity is the only force acting at a distance on it.

Boulton ([c1975]) derives a condition for the lodgement of a stone on a glacier bed by estimating the horizontal force exerted on the stone by the ice and equating this to the retarding friction between stone and bed. Strictly speaking, the analysis in this paper cannot be applied to a stone in contact with the bed. However, if the stone is extremely close to the bed, i.e. when ∆ ≈ R, we might postulate the existence of a very small irregularity on the stone which allows normal and frictional forces to be transmitted without disturbing the pattern of flow of the ice or altering the physical processes at the boundaries. Then we might hope to derive a more accurate estimate of the lodgement condition, based on an analytical solution for the forces of the ice on the stone rather than an order-of-magnitude calculation.

The Navier–Stokes equations in bipolar coordinates

The bipolar coordinates α, β are defined in terms of the Cartesian coordinates x, y by the equations


The coordinate curves are two sets of coaxial circles (Fig. 2) which intersect orthogonally. The scale factor a is the distance between the point α = ∞ and the line α = 0 (the x-axis). The metric coefficient is


Elements of arc have lengths δα/h, δβ/h. The radius of the circle α = α 1, is R = a/sinhα 1, and the distance of its centre from the line α = 0 is ∆ = a|coth α 1|. In this coordinate system, the Navier–Stokes equations for the slow steady flow of a linearly viscous incompressible fluid under gravity are




Fα and Fβ are the components of the body force which can be derived from the hydrostatic pressure p 1


If the line α = 0 is horizontal, then


where p 1 = p b is constant on α = 0, ρ is the density of the ice, and g the acceleration due to gravity.

Fig. 2. The bipolar coordinate system.

The components of the stress tensor σ ij are related to the components of the strain-rate tensor ė ij by the equations


where η is the viscosity and p is a pressure term. The continuity equation is


where uα is the velocity in the α direction and uβ the velocity in the β direction. Thus, a stream function ψ may be defined by the equations


The Navier–Stokes and continuity equations reduce to the biharmonic equation for the stream function


and the Laplace equation for the pressure p


The general solution for the biharmonic equation in bipolar coordinates has been given by Jeffrey (1920, 1922)


where the Fourier coefficients ɸn , Χn are functions of α:


Since p and η2 Ѱ are conjugate functions, an expression for p in terms of the constants etc., can be derived. The leading term is 2η(B 0+D 1) β/a. Since p must on physical grounds be single-valued with β,


The components of the strain-rate tensor can be expressed in terms of ψ


and hence, using Equations (12) and (13), may also be written as Fourier series in β with coefficients that are functions of the coefficients An, An ′, etc. From Equation (6), p 1 may also be expanded in a Fourier series


Hence, from Equation (7), the components of the stress tensor may be expanded in Fourier series with coefficients that are functions of the coefficients.

Boundary conditions

Figure 1 shows a cylindrical stone, radius R, at a distance ∆ from the horizontal bed of a glacier. The surfaces of the stone and bed are defined by the curves α = α 1, α = 0 respectively. The equation of the upper surface of the ice is


Up- and down-stream boundaries are defined by the equation


The coefficients in Equations (13) are determined by applying boundary conditions to these five curves.

(i) The ice–rock interfaces

Conditions at these interfaces will be described in terms of the “classical” regelation theory (Nye, 1967) although laboratory experiments have shown that this theory is oversimplified (Nunn and Rowell, 1967; Townsend and Vickery, 1967; Drake and Shreve, 1973; Morris, 1976). Drake and Shreve have pointed out that the effect of impurities, supercooling, and the formation of a trace must be included in a complete regelation analysis. However, we use the simple theory since it allows a simple and clear mathematical demonstration of a difficulty which I believe would arise in any case. In this geometry the classical regelation analysis proceeds as follows: Between ice and rock on the surfaces α = o, α = α 1 there are thin films of water, the regelation layers. On the boundary between the ice and the stone (α = α 1) there is no shear stress;


On part of the boundary between the ice and the bed we may want to allow a non-zero shear stress τ produced by the effective roughness. Thus


where β 1 and β 2 are arbitrary constants.

This boundary condition can be written in the form of a Fourier expansion for σαβ,


on the surface α = 0.

In areas where the normal pressure on the rock is high, ice melts and the water travels through the regelation layer to the areas of low normal pressure. The ice–water interfaces are always at the melting-point temperature, which varies with the stress normal to the interface. Thus


where C is a constant having the value 0.7 x 10–7 deg Pa–1.

The high-pressure areas have a lower temperature than the low-pressure areas. Thus the latent heat released by the ice when it refreezes in the low-pressure areas flows through ice, stone, and bed to melt more ice in the high-pressure areas.

The velocity of the ice in the positive direction is


The latent heat released by freezing per unit time on an element δβ/h is Luαδβ/h on α = 0 and —Luαδβ/h on α = α 1, L is the latent heat of melting of ice per unit volume (2.8 X 108 J m–3). Let T i, T s, T b be the temperature distributions in the ice, stone, and bed respectively, and k i, k r the thermal conductivities of the ice and rock. Then, on α = 0,


and on α = α 1


In the classical regelation analysis the temperature distributions are solutions of the Poisson equation with boundary conditions given by Equations (24) and (25), and the restriction that the temperature should remain finite as α → 0, β → 0 (i.e. as x → ∞, y → ∞) and as α → ± ∞.

The distributions have the form




Using Equation (22) the T coefficients, En , Gn , etc., can be denned in terms of the coefficients. Then Equations (24) and (25) give sets of equations relating groups of the coefficients. Further equations are derived from the boundary conditions (19) and (21) for the shear stress. We arrive finally at sets of linear equations from which all the coefficients may be determined given A 0, B 0, C 0, D 0, B 1, A 2, C 2, A 3, C 3, A 1′, B 1′, C 1′, D 1′, A 2′, and C 2′. The further equations which determine the values of these last coefficients are derived from other boundary conditions.

(ii) The ice-air. interface

Since Ha we have α → 0, β→ 0 on the upper surface of the glacier. The expansions for p and ėαα reduce to

Therefore, from Equation (7) the boundary condition

reduces to


Now h → 0 as α and β → 0 so huα 0, → 0 if uα and ѱ remain finite on y = H. We may therefore specify two more conditions


(iii) The up-stream and down-stream boundaries

In a complete analysis of glacier flow the gravitational forces which produce motion in the ice are determined from the slope of the bed and/or the upper surface of the ice (e.g. Morland, 1976). We have focussed attention on a small horizontal area of the bed over which the ice has a uniform thickness. We have to suppose that up-stream and down-stream of this area gravitational forces act so as to produce a certain pattern of flow of the ice as it enters and leaves the area of interest around the stone. This pattern is defined by boundary conditions on x = ±X 1. For example, we might suppose that far away from the stone there is simple shear flow parallel to the bed. Thus, for example,


where u b is the sliding velocity at the bed.

Forces on the stone

Let the forces applied to the surface α = α 1 be statically equivalent to forces X and Y acting at the centre of the stone and a couple L. X acts in the direction parallel to the bed and Y in the perpendicular direction


On the surface of the stone σα β = 0. From Equations (22) and (26), using the series expansions


Expressions for the forces may be obtained in terms of the T coefficients which are themselves known in terms of the coefficients. Thus,



The solution of the equations relating the coefficients of the series expansion for (Equation (12)) is discussed in the Appendix. It is found that the series of coefficients {ɸn } and {ѱn } diverge, unless both λ and λ′ (two dimensionless critical wavelengths) are zero. Given that L, C, k r, k i, and η are known, finite, and non-zero physical constants, this condition means that a → ∞ and the ratio of the distance of the radius of the stone to the distance of its centre from the bed R/ ∆ → 0. In other words, a solution exists for an isolated cylinder in ice (and has indeed been given by Watts (unpublished)), but not when there is another solid boundary.

Thus, we have found that, in the steady state, the classic regelation boundary conditions cannot be obeyed both at a flat bed and at the surface of a circular stone. That is, the temperature and stress distributions cannot be matched so that both the surfaces are at the pressure-melting point. A similar problem has arisen before in theoretical work on regelation and has been discussed by Nye (1967) and Morris (1976). In their analyses, for wires and spheres moving through ice and for ice moving around a cylinder with a wavy surface respectively, there is apparently only one boundary at which normal stress and temperature must be matched. However, at any point within the ice where there is a water inclusion, for example at a three-grain intersection, a relation between, the normal stress across the interface and the temperature can also be defined. Thus there are internal boundaries to be considered. Morris suggested that melting and refreezing within the ice would produce an “internal” temperature distribution to be added to the “regelation” temperature distribution which is produced by melting and refreezing on the solid boundary. The total temperature distribution would match the stress distribution so that the ice is always at the pressure-melting point. The classic method of analysis of regelation problems depends on the assumption that any “internal” component of temperature is negligible compared to the “regelation” component at the solid boundary.

I suggest that, even in the time-dependent case, when there are two solid boundaries with melting and refreezing on each there is no reason to suppose that the temperature and stress distributions can be matched so that the ice at both surfaces is at the pressure melting point. There are three ways in which this problem could be resolved:

  • (1) melting and refreezing within the ice could produce an “internal” temperature distribution such that the ice at the two boundaries and at the internal water inclusions was at the pressure-melting point. The magnitude of the internal heat sources and sinks would not be negligible compared to the magnitude of the sources and sinks on the boundaries;

  • (2) the ice could separate from the boundaries producing water-filled cavities with shapes such that the temperature at the ice–water interface was always the correct melting-point temperature for the normal stress across the interface. At the moment we do not know if there is any shape for which a steady-state solution is mathematically possible given our formulation of the problem with linearly viscous ice and classical regelation boundary conditions;

  • (3) the ice could separate from the boundaries producing cavities filled with air and/or water vapour. In this case, parts of the boundary where there is no melting can be well below the pressure-melting point and the requirement that stress and temperature distributions should match at all points of the boundaries is relaxed.

Field observations of clasts in basal ice (e.g. Vivian and Bocquet, 1973; Boulton and others, 1979) indicate that cavities frequently occur in the lee of large particles. Some of these cavities are full of water, others have long spicules of clear ice. The habit of these ice crystals suggests that they may have grown into a cold atmosphere which is below the pressure-melting point.

Of course, the conditions at the bed of a glacier are vastly more complicated than the simple problem that has been analysed in this paper. We have, for example, ignored the effect of geothermal heating, drainage of water at the glacier bed, and the interaction between several clasts in transport which may obscure the effects of the stress–temperature adjustment process described here. The classical regelation theory used in this paper and in current basal sliding and lodgement theories is known to be inadequate. However, two conclusions can be drawn which are important for theories of erosion and deposition:

  • 1. Cavities may form around clasts in transport in the basal ice because of the presence of the bed. These cavities must be distinguished from those produced by “cavitation” (Nye, 1970) which is controlled by the overburden pressure of ice.

  • 2. The size and shape of such cavities will depend on the pattern of flow of the ice. If the stone moves to a part of the bed with a different flow pattern, or the pattern at a given point changes in time, bulk melting or refreezing must take place as the cavities adjust to the new situation.

Thus, we cannot expect that the forces on a stone due to the ice can be calculated by a simple application of classic regelation and plastic-flow theory. Even when the clast is in contact with the bed and there is only one ice–rock boundary there may not be solutions of the partial differential equations for ѱ and T if regelation boundary conditions are applied along the whole of the boundary. Happel and Brenner ([c1973], p. 61) remark that it is difficult to generalize on the required conditions for the existence of unique solutions to the Navier–Stokes equations given combinations of prescribed velocities and derivatives at the boundary. We know there are solutions for the rather odd mixed boundary conditions that arise from classical regelation theory when these are applied on spherical, cylindrical, and plane surfaces (note that in Nye’s solution for a perturbed plane and the Morris adaptation of this for a perturbed circular cylinder the boundary conditions are applied at the unperturbed surface). However, these are rather special geometries in which the orthogonal curvilinear coordinate system defined by the boundary curve has metric coefficients which do not vary with position along the boundary. Whether solutions exist for Newtonian flow of ice over rock humps or clasts of any shape resting on the bed with classic regelation conditions applied over the whole boundary is an open question.

I do not believe that the problem of the adjustment of stress and temperature distributions that has been raised in this paper is an artefact of the particular geometry, flow law, and boundary conditions chosen. The simplifications I have made merely allow a clear analytical demonstration of the need for an adjustment process. This process will almost certainly occur whenever there are two ice–rock interfaces close together and may also have a part to play in a complete description of flow over steep bedrock obstacles.


This Appendix gives the derivation of the coefficients of ɸ n for high n. Similar equations hold for the coefficients of Xn . For low values of n some terms of the equations are lost or altered but the overall form is the same.

The shear-stress boundary conditions (Equations (21) and (19)) lead to




The normal stress boundary conditions lead to equations for the T coefficients


Ri and Si are functions of n, cosh 1 and sinh 1:


If ux /h expanded in the Fourier series




the velocity boundary conditions (Equations (24) and (25)) give




Substitution for Fn and Hn from Equations (A-3) and (A-4) into Equations (A-7) and (A-8) gives




λ and λ′ are dimensionless critical wavelengths which may be compared to the critical wavelength given by Nye (1969):




Now, from Equation (9)


hence from Equations (A-5) and (A-6)




Substitution for ρ i and ϒ i from Equations (A-9) and (A-10) and for B i and D i from Equations (A-1) and (A-2) using i = n — 3, n+3 leads to two simultaneous equations for A n+3 and C n+3 in terms of the known coefficients A n-3 to A n+2 and C n-3 to C n+2. Note that some coefficients are functions of λ or λ′, others depend only on α 1.


Boulton, G. S. [c 1975.] Processes and patterns of subglacial sedimentation: a theoretical approach. (In Wright, A. E. and Moseley, F. ed. Ice ages: ancient and modern. Liverpool, Seel House Press, p. 742. (Geological Journal Special Issue No. 6.))
Boulton, G. S. and others. 1979. Direct measurement of stress at the base of a glacier, by Boulton, G. S. Morris, E. M. Armstrong, A. A. and Thomas, A.. Journal of Glaċiology, Vol. 22, No. 86, p. 324.
Drake, L. D. and Shreve, R. L. 1973. Pressure melting and regelation of ice by round wires. Proceedings of the Royal Society of London, Ser. A, Vol. 332, No. 1588, p. 5183.
Glen, J. W. and others. 1957. On the mechanism by which stones in till become oriented, by Glen, J. W. Donner, J. J. and West, R. G.. American Journal of Science, Vol. 255, No. 3, p. 194205.
Happel, J. and Brenner, H. [c 1973.] Low Reynolds number hydrodynamics with special applications to particulate media. Second revised edition. Leyden, Noordhoff International Publishing.
Jeffery, G. B. 1920. Plane stress and plane strain in bipolar co-ordinates. Philosophical Transactions of the Royal Society of London, Ser. A, Vol. 221, No. 9, p. 26593.
Jeffery, G. B. 1922. The rotation of two circular cylinders in a viscous fluid. Proceedings of the Royal Society of London, Ser. A, Vol. 101, No. 709, p. 16974.
Kamb, W. B. 1970. Sliding motion of glaciers: theory and observation. Reviews of Geophysics and Space Physics, Vol. 8, No. 4, p. 673728.
Lewis, W. V. ed. 1960. Investigations on Norwegian cirque glaciers. London, Royal Geographical Society. (R.G.S. Research Series, No. 4.)
Lliboutry, L. A. 1968. General theory of subglacial cavitation and sliding of temperate glaciers. Journal of Glaciology, Vol. 7, No. 49, p. 2158.
Morland, L. W. 1976. Glacier sliding down an inclined wavy bed. Journal of Glaciology, Vol. 17, No. 77, p. 44762.10.1017/S0022143000013733
Morris, E. M. 1976. An experimental study of the motion of ice past obstacles by the process of regelation. Journal of Glaciology, Vol. 17, No. 75, p. 7998.
Nunn, K. R. and Rowell, D. M. 1967. Regelation experiments with wires. Philosophical Magazine, Eighth Ser., Vol. 16, No. 144, p. 128183.
Nye, J. F. 1967. Theory of regelation. Philosophical Magazine, Eighth Ser., Vol. 16, No. 144, p. 124966.
Nye, J. F. 1969. A calculation on the sliding of ice over a wavy surface using a Newtonian viscous approximation. Proceedings of the Royal Society of London, Ser. A, Vol. 311, No. 1506, p. 44567.
Nye, J. F. 1970. Glacier sliding without cavitation in a linear viscous approximation. Proceedings of the Royal Society of London, Ser. A, Vol. 315, No. 1522, p. 381403.
Röthlisberger, H. 1968. Erosive processes which are likely to accentuate or reduce the bottom relief of valley glaciers. Union de Géodésie et Géophysique Internationale. Association Internationale d’Hydrologie Scientifique. Assemblée générate de Berne, 25 sept.–7 oct. 1967. [Commission de Neiges et Glaces.] Rapports et discussions, p. 8797. (Publication No. 79 de l’Association Internationale d’Hydrologie Scientifique.)
Townsend, D. W. and Vickery, R. P. 1967. An experiment in regelation. Philosophical Magazine, Eighth Ser., Vol. 16, No. 144, p. 127580.
Vivian, R. A. and Bocquet, G. 1973. Subglacial cavitation phenomena under the Glacier d’Argentière, Mont Blanc, France. Journal of Glaciology, Vol. 12, No. 66, p. 43951.
Watts, P. A. Unpublished. Inclusions in ice. [Ph.D. thesis, University of Bristol, 1974.]
Weertman, J. 1957. On the sliding of glaciers. Journal of Glaciology, Vol. 3, No. 21, p. 3338.


J. Weertman: Do your boundary conditions take into account the surplus, unfrozen water that is left over after the regelation cycle is completed? The volume of this water produced per unit time should be equal to the work done per unit time divided by the latent heat of melting–freezing per unit volume of ice. You could not obtain a consistent set of equations that could give a steady-state solution for this problem. But if water is continuously produced by the regelation cycle—which is a Carnot cycle—the problem cannot be a steady-state one.

E. M. Morris: My boundary conditions are those used by other authors such as G. S. Boulton and B. Hallet, since I am trying to investigate the validity of their order-of-magnitude analysis. I agree that the regelation equations used by these authors are inadequate.