Ice Flow Over a Bump: Experiment and Numerical Simulations

Abstract

The flow of ice over a bump under steady-state conditions has been studied both experimentally and numerically. The experiments were performed with the ice viscometer installed at the Laboratoire de Glaciologie et Géophysique de l'Environnement du C.N.R.S. (LGGE) in Grenoble, France. The flow was observed by means of a marker network photographed at fixed time intervals. The separation of the ice down-stream of the bump and the non-linearity of the generally accepted Glen’s flow law led to the modelling of this flow using a finite-element method. The results describe the velocity and stress fields in ice moving over a bump and obeying Glen’s flow law (until now only relatively rough approximations of these fields have been given in the literature). The comparison of the observed and calculated network deformations shows that, although ice flowing over a bump does not undergo steady creep, the use of Glen’s flow law with an exponent n = 3 remains justified when simulating the flow of a glacier over its bed under steady-state conditions.

1. Introduction

Research on glacier sliding laws has been the subject of considerable work and numerous controversies since Weertman’s first model (Weertman, 1957, 1964, 1971; Lliboutry, 1968, 1975, 1979; Nye, 1969; Kamb, 1970). Even if a very simple micro-scale model is used for the glacier bed, several major difficulties remain in describing the flow of ice over an isolated obstacle. For the small obstacles on the glacier bed, two main phenomena are involved: (1) sliding by “melting-refreezing”, the simulation of which still poses theoretical and experimental problems as the type of surface in contact with the ice, its salt and air-bubble contents, its permeability, and the way in which the liquid-water inclusions migrate must be accounted for (Drake and Shreve, 1973; Nye, 1973; Chadbourne and others, 1975; Lliboutry, 1976; Morris, 1976; Shreve, 1984); (2) the deformation of ice, difficult to integrate in a sliding law, mainly because of the non-linearity of the generally adopted constitutive law. These two processes are often accompanied by separation of the ice down-stream of the obstacles (“cavitation”) which complicates the problem even further because: (1) the flow over an isolated bump can only be calculated numerically, even if a Newtonian viscous law is adopted for ice (Iken, 1981); (2) to account fully for cavitation in a sliding law, a model giving the subglacial water pressure is required (Röthlisberger, 1972; Nye, 1976; Iken and others, 1983; Lliboutry, 1983).

The present work, based on experiments carried out with the “Penelope glacier simulator” located at LGGE-Grenoble, is only a small contribution to the problems of establishing a glacier-sliding theory and lies in the wake of studies by Brepson (1979, 1980) and Hooke and Iverson (1985). It consists of a localized study of the flow of ice over a succession of identical bumps, showing a strong separation down-stream of these bumps. The melting-refreezing process is ignored and the down-stream cavity water pressure is assumed to be known. We restrict ourselves here to checking whether Glen’s flow law, which is generally accepted to describe the steady creep of ice, can be used to simulate the flow of a glacier near its bed under steady-state conditions. This requires the numerical simulation of an experiment selected for its long duration.

2. Description of the Experiment

2.1. Apparatus

A detailed description of the “Penelope glacier simulator” has been given by Lliboutry and Brepson (1963) and Brepson (1979). Only the essential features of the apparatus, as well as the modifications made to the original device, are described here.

The ice flow takes place in a large-diameter Couette-type viscometer (see Fig. 1). The outer cylinder (1) is a notched ring, made of an insulating material, which sets the ice in motion. The inner cylinder (2), made of epoxy resin, has an almost elliptical cross-section and represents two identical smooth bumps placed symmetrically about the minor axis. The sliding velocity of the ice relative to the bumps can be set from 1 to 1000 m a⁻¹. The inner cylinder (2) is connected to the fixed frame of the machine by means of a calibrated torsion rod, the rotation of which is monitored using a LVTD. The lower part of the flow chamber is delimited by an insulating plate (3) which turns with the driving ring (1). The upper part of the apparatus, described by Brepson (1979), has been eliminated and replaced by a thick Plexiglass plate (4) which showed negligible deformation during the tests. The outer part of this plate is fastened to the notched ring (1); its inner edge rests against the top of the inner cylinder (2) via a flat sliding pad. Contrary to the original design, the upper plate (4) cannot move in a vertical direction. The ice sample is cog-wheel shaped and has an outer diameter of 60 cm, a mean inner diameter of 37 cm and a thickness of 13 cm.

Fig. 1. Schematic diagram of the apparatus; (1) notched driving ring, (2) bumps, (3) lower plate, (4) upper transparent plate.

The two horizontal plates (3) and (4), which grip the ice and are rotated with it, have smooth surfaces and, since they are lubricated by a thin water film, a frictionless contact with the deforming ice can be assumed. Under these conditions, plane flow is achieved.

Finally, a built-in cooling circuit allows the temperature of the ice to be modified without changing the cold-room setting.

2.2. Ice-sample preparation

The ice was made from a mixture of snow and water. Snow was put in the flow chamber then soaked in deionized water which was kept at 0°C This mixture, which contained about 50% snow and 50% water by weight, was stirred manually then refrigerated for several hours using the internal cooling circuit. Ice obtained in this way has an isotropic structure (Duval and Ie Gac, 1980). Its density, calculated from the masses of snow and water used, was about 0.83. The cooling system was then stopped and, when the sample reached cold-room temperature, the machine was set in rotation.

During the experiment described here, ice densification was such that the bumps were completely cleared after one revolution of the driving ring. The resulting cavities were filled up with snow and water and the process was then repeated. At the end of the experiment, the exact density was found to be 0.89 ± 0.01 (average on four ice pieces taken from the sample).

2.3. Flow visualization

The flow was observed using an array of uniformly distributed markers (see Fig. 2). The markers were polyethylene rods, 3 mm in diameter and 7 cm long, planted vertically in the ice. The ice pressure exerted on their lower section flattened their upper section against the Plexiglass plate, making them easy to distinguish throughout the whole experiment. A camera took photographs of the network at regular intervals of 63 min. It was also possible to follow the evolution of the cavities formed down-stream of the bumps, from the beginning of the experiment (see Fig. 3).

2.4. Experimental procedure

As said previously, physical modelling of the melting-refreezing process remains difficult. On the other hand, even using the classical simplified theory (i.e. heat transfer by pure conduction through rock and ice; see, for example, Lliboutry, 1975), the relatively small thickness of the insulating material forming the bumps and the horizontal plates (2—5 cm), and the surrounding metal parts, would lead to a three-dimensional coupled problem. The solution of such a problem lies outside the scope of this study and thus the influence of the melting-refreezing process must be limited, if not eliminated. Consequently, the experiment was carried out at a temperature between −0.5° and −1°C: the ice must be slightly below the melting point corresponding to the mean normal stress acting at the ice-bump interface estimated from the torque measurements.

During the first few days the speed of rotation was set very low to prevent the torque from reaching the torsion-rod safety limit (C _max = 5500 Nm). This period corresponds to the opening phase of the cavities by ice densification. Since the expansion of the cavities resulted in a decrease of the torque, the speed of rotation was increased progressively so that the torque reached 0.9C _max and then set to a fixed value of around 1/4 revolution/d until the end of the experiment. When the ice density stopped increasing and the calorific balance was achieved, the down-stream cavities stabilized and the ice-bump contact area remained almost constant. Nevertheless, the torque continued to oscillate slowly around its mean value C _m, with a frequency of about 2 cycles/24 h and an amplitude of 0.05C _m (see Fig. 4). Once this virtually steady state was reached, the experiment was continued for as long as possible to obtain a sufficient deformation of the marker array. In all, the experiment lasted 29 d, of which only the last 16 d showed acceptable steady-state conditions. Time t₀, in Figure 4, which will be referred to as the “initial time”, corresponds to the beginning of this quasi-steady-state period (and to the end of the first revolution at fixed speed of rotation).

Fig. 4. Measured torque versus time. (Steady state was assumed at t = t₀; one revolution between t₀ and t₁ four revolutions between t₀ and t₄.)

3. Numerical Simulation

3.1. Description of the problem

3.1.1. Equations

The melting-refreezing process is ignored. The ice is assumed to be a continuous, homogeneous, isotropic, incompressible, isothermal medium. Disregarding the effects of gravity and inertia, the equations of equilibrium and mass conservation are:

(1)

(2)

where σ_ij are the stress tensor components and u is the velocity vector.

3.1.2. Constitutive law for ice

The assumed constitutive law is Norton-Hoffs power law which relates the strain-rates έ _ij to the deviatoric stresses σ’_ij:

(3)

where τ is the effective shear stress defined by: This constitutive law leads to the relation known as Glen’s flow law (Glen, 1955): where is the effective shear strain-rate given by

The steady-creep law, equation (3), was verified by Duval (1976) in torsion-compression tests on natural ices (n was found equal to 3 for τ in the range 0.1–0.4 MPa). It was used in our computations, although when passing over a bump a material particle does not undergo steady-creep conditions (Paterson, 1981). The isothermal ice hypothesis led to the use of a constant for B.

3.1.3. Coordinate systems, flow geometry, notation

The numerical solution is achieved using a Cartesian reference frame with the x-axis joining the summits S-S’ of the bumps and the z-axis along the rotation axis of the machine (see Figs 1 and 5). Plane flow is assumed in the x—y plane. For easier reading, the results are shown in polar coordinates (ρ,θ) as defined in Figure 5.

Fig. 5. Mesh used for computation, coordinates, and notation (shadowed area: ice-bump contact)

The geometry of the flow domain D is simplified by replacing the actual notched outer boundary with a circle of equation ρ = R_e (see Fig. 5). As steady flow is assumed, the inner boundary formed by the two cavity roofs and ice-bump interfaces is a stream line and a trajectory with an equation ρ = ρ₀(θ). In the following, the point at which an ice particle moving along this trajectory loses contact with a bump is named the “separation point”, the point at which it re-contacts the next bump is named the “rejoining point”.

The following notation and conventions are used:

u Velocity vector.

u_x, u_y Velocity components in rectangular coordinates.

u _ρ , u _θ Velocity components in polar coordinates.

n Unit normal to the boundary of D, pointing outward from the ice.

t Unity tangent obtained by anti-clockwise rotation of n by π/2.

Σ Stress vector acting on the boundary of D.

Σ_n, Σ_t Normal and tangential components of Σ, along n and t (positive when tensile)

p Isotropic pressure (positive when tensile).

Ω Angular (driving) velocity of the apparatus.

3.1.4. Boundary conditions

Taking into account the symmetry with respect to the z-axis, the solution, expressed in polar coordinates, is periodical of period π. Consequently, the modelled region is halved (see Fig. 5).

The prescribed boundary conditions are as follows:

• (a) periodicity π expressed by setting u(x,0) = −u(−x,0) and Σ(х,0) = Σ(−x,0) (on ST and S’ T’ boundaries in Figure 5);

• (b) constant angular velocity Ω, resulting in tangential

  velocity U of constant magnitude on the outer ring TT’ : u_ρ = 0; u _θ = U = R _eΩ

• (c) free-surface condition Σ = 0 on the cavity roof S’ C;

• (d) frictionless ice-bump contact expressed by: un = 0 and Σ_t = 0 on SC. It must be verified a posteriori that the ice-bump interface is subject only to compression (Σ_n > 0).

3.1.5. Comments

Condition (c) above corresponds to a zero air pressure P_c inside the cavity. The solution for p_c = H = constant is simply obtained by adding H to the isotropic pressure (the velocity and deviatoric stress fields are unchanged). As ice incompressibility is assumed, the volume of the cavity, fixed by the initial shape given to the domain D, does not depend on p_c. The problem is thus different from the one treated by Brepson (1979, 1980): in his experiments the upper plate was free to move vertically and its upper side was submitted to water pressure p_w. During steady plane flow, the shape of the cavity adapted itself so that the action of the vertical stresses exerted by ice balanced the action of p_w.

3.1.6. Choice of units

The computation is made using dimensionless variables. The length unit is denoted by L. The magnitude U of the tangential driving velocity is taken as the velocity unit. Consequently, the strain-rates are expressed in units of U/L. For the sake of convenience, the Glen’s flow law coefficient B is a fixed a priori. The stress unit is therefore σ = (U/BL)^1/n .

As the reduced equations of the problem are identical to Equations (1) and (2), the notation used for the actual and reduced variables is the same.

3.2. Computation method

3.2.1. Variational formulation

The problem is put into a form suitable for numerical solution by using the Bird (1960) Johnson (1960) variational principle for steady flow of a non-Newtonian fluid with dissipation potential Φ. It has been proved (Meyssonnier, 1983) that the solution for a given cavity shape is obtained by making the functional

stationary over the set of kinematically admissible velocity fields, i.e. those which satisfy here:

• (a) u(x, 0) = −u(−x, 0);

• (b) u = U on the outer ring (ρ = R _e, 0 ≼ θ ≼ π);

  where uⁱ _j, pⁱ are the values of u_x, u_y, and p at the Ne and Me nodes i and are polynomial functions of x and y.

  The stationarity of J is expressed by solving the set of simultaneous equations with the unknowns uⁱ _j, pⁱ

• (c) un = 0 on the ice-bump contact.

3.2.2. Finite-element method solution

This formulation leads naturally to the choice of velocity and pressure as primitive variables (Nickell and others, 1974; Thompson, 1975). The solution to δJ = 0 is carried out by discretizing domain D into a finite number of geometrically simple elements. In each element e the velocity U^e and pressure p^e are interpolated as:

(4)

where N is the total number of “velocity nodes” and M is the total number of “pressure nodes”, while prescribing boundary conditions on the velocities.

In this work Lagrangian triangular six-node elements, with a quadratic interpolation of the velocity and a linear interpolation of the pressure, were chosen (Thompson, 1975; Jackson and Cliffe, 1981). Continuity of the velocity across triangle sides is ensured, but not that of the pressure so as to force complete incompressibility at each point within each element. So each triangle has six “velocity nodes” (three vertices + three mid-side nodes) and three internal nodes corresponding to the “pressure at triangle vertices” unknowns. The mesh was made up of 142 elements and 351 “velocity” nodes, of which 105 were triangle vertices (see Fig. 5). Its regular structure was broken so as to refine the element size near the rejoining point. Curved triangles were used to describe the inner and outer boundaries of D (strict incompressibility in these elements is then no longer ensured). It is worth pointing out that this mesh discretizes the space occupied by ice and not the ice itself, corresponding to an Eulerian description of the flow.

3.2.3. Assembly and solution of the finite-element system

Equations (4) are assembled using the contributions from each element to integral J. The boundary conditions on u are specified after assembly. The driving velocity U is prescribed using the Payne—Irons method. The sliding and periodicity conditions are prescribed by making the appropriate combinations of rows and columns of the system, and eliminating the excess unknowns. The resulting non-linear finite-element system is:

(5)

where Y is the nodal velocity and pressure-unknown vector, the matrix K depends on the velocities through the viscosity and B ₀ is the “load” vector resulting from the prescribed velocity condition. It is solved by successive updating, for which the convergence, for n ≥ 1 in Glen’s flow law, was proved by Friâa (1979). At step i the viscosity at each integrating point is calculated using the strain-rates obtained at step i − 1 as The resulting linear system, [K(u _{i − 1})]Y_i = B ₀, is solved by the Gauss algorithm.

In this work, the iterations were stopped at step i when ||Y_i − Y _{i − 1}|| / ||Y_i < 10⁻³ and ||B ₀ − [K(u_j )]Y_i ||² < 10^–5 (σL ²)², where ||V||² is the squared length of vector V.

3.2.4. Free-surface fitting

At each step the solution verifies only the natural condition Σ = 0 on the cavity roof S’ C (see Fig. 5): the extra requirement for steady flow is that S’ C must be a stream line. This is carried out as follows:

• (a) the initial shape of the cavity is fixed a priori as the arc of circle tangent to one bump at its summit S’ (ρ₀, π) and intersecting the other bump at point C (ρ₀, π/4), close to the observed position;

• (b) after calculating the solution to system (5), the trajectory of the ice particles passing through S’ is obtained by integrating the equation dy/dx = u_y/u_x;

• (c) this trajectory is then taken as the new free surface and the corresponding velocity field is re-computed;

• (d) the different computations are re-iterated until the condition un = 0 is satisfied on the free surface.

At each stage one must verify that the normal stresses at the ice-bump interface nodes are all compressive, and that the new computed free surface does not overlap the bump profile. This represents the main difficulty of the process, due to the changes in nature of the boundary conditions occurring near the ends of the cavity.

From the start, it appeared that the ice tended to contact with the bump down-stream, looking like the swelling of jets in extrusion die outlets (Nickell and others, 1974). The iterations were stopped when the ratio |un/uu| was lower than 10⁻² at all the nodes on the free surface. Their displacements, normal to the cavity roof, were then below 10⁻⁴ L.

3.3. Estimation of strain-rates and stresses

The interpolation used for the velocities results in strain-rates έ _ij which are discontinuous across the triangle boundaries. In order to obtain easily usable information, a continuous description, in terms of nodal values, of the quantities presenting such discontinuities is needed. Many smoothing techniques are available (see, for example, Lee and others, 1979). The simplest method, consisting of calculating local nodal averages, was used in this study: being the value calculated at node i in triangle e, and N the number of triangles including node i, the smoothed nodal value is simply

The mean pressures at the triangle vertices were calculated in this way, using the elementary values obtained as part of the solution to system (5). Then the pressure field was interpolated linearly through the vertices. The elementary values of the strain-rates, viscosity, and deviatoric stresses were calculated at the six nodes of each triangle, then smoothed to obtain nodal values. Suitable descriptions of the respective fields were obtained by adopting the quadratic interpolation used for the velocities. This has been assessed in the case of the flow of a Newtonian material, without cavitation, for which an approximate analytical solution may be obtained (Meyssonnier, 1983). In the present case, it is only possible to check the validity of the results:

• (a) on the inner stream line, the computed stresses correctly account for the prescribed boundary conditions, except in the immediate vicinity of the separation and rejoining points, which present singularities (see Fig. 6);

• (b) on the outer ring (ρ = R_e ), the condition u _ρ = 0; u _θ = U= constant leads in theory to έ_ρρ = 0 (Meyssonnier, 1983): this condition is achieved with , except at the ends (θ = 0,π) where ;

• (c) the value of the dissipated power computed from the strain-rates, , and that derived from the torque C exerted by the normal stresses acting on the bump P_e = CΩ are consistent since: (P _d − P_e )/P_e = 3 *#x00D7; 10^-2.

Fig. 6. Normal (Σ_n) and tangential (Σ_t) components of stress acting on the ice at the bump contact and on the cavity roof (dotted line: isotropic pressure p)

4. Numerical and Experimental Results

4.1. Numerical results

4.1.1. Velocity field

Due to the large extent of the cavity, the velocity field is little disturbed in the zone situated above it. The u _θ profiles are close to those of a rigid body in the sector π/2 > θ > π(see Fig. 7). The disturbance becomes important in the vicinity of the rejoining point.

Fig. 7. Profiles of the u_θ tangential velocity component along radii of polar angles θ = 0, θ = π/4, θ = π/2, computed with n = 3 (dotted lines: rigid motion u_θ = ρΩ profiles).

4.1.2. Strain-rates

Contours of derived from the computed velocity field are shown in Figure 8. This figure can be read as a map of τ contours by applying Glen’s flow law. The strain-rate gradient is very high in the region near the rejoining point where γ passes from 2 to 17.6 (U/L)its maximum value, over a distance equal to 0.015L, giving a gradient of about 10³(U/L ²)•

Fig. 8. Effective shear strain-rate contours (Glen’s flow law, n = 3; dimensionless values).

4.1.3. Stresses

The nodal isotropic pressures are all compressive, lying between −0.2σ and −7.7σ. The minimum compression is reached on the cavity roof, just after separation. The maximum compression is reached at the rejoining point. The istropic pressure remains high in the whole region of convergence of the flow (see Fig. 9), being very high at the rejoining point and stretching over the ice-bump interface. Thus, the isotropic pressure contributes largely to the normal stress acting on the bump (see Fig. 6). A similar result was found by Lliboutry and Ritz (1978) for ice flow past a sphere. Figure 10 shows the principal stresses computed at the centroids of the triangles. All of these stresses are compressive. The influence of the free surface is significant and remains perceptible in the core of the ice above the cavity. On the other hand, the extent of the cavity-bump transition zones within the ice, characterized by the rotation of the principal stress directions, seems very limited.

Fig. 9. Isotropic pressure contours (Glen’s flow law, n = 3; dimensionless values).

Fig. 10. Principal stresses at the element centroids (n = 3).

4.2. Relation between computation and experiment

4.2.1. Units and value of B in Glen’s flow law (3)

The length unit such as the dimensions of the finite-element mesh correspond to that of the apparatus is L = 1 m. The angular velocity O was measured very accurately thanks to the photographs and checked to remain constant for the duration of the experiment. It corresponded to a driving velocity U = 154.1 ± 0.5 m a⁻¹. The resulting strain-rate unit is U/L = 154 a⁻¹. The stress unit is derived from the comparison between the measured and computed values of the torque exerted by the ice on the torsion rod. The “theoretical” torque calculated from the Σ_n distribution shown in Figure 6, with a sample thickness of 0.13L is C _th = 48.5 × 10⁻⁴ (σL³). The “experimental” torque is taken as the average of the observed values (see Fig. 4), T being the period during which quasi-steady state was achieved. According to Mellor and Cole (1981), a better value C _e, accounting for the non-linearity of Glen’s flow law, would be such that

In the present case, the fluctuations of the torque, compared to C _e or C _exp, are small enough so that the two estimates are very close. The average value, taking into account a friction torque of 600 N m, is C _exp = 4150 Nm The C _exp/C _th ratio gives the stress unit: σ = 0.85 Mpa (8.5 bar).

The resulting viscosity unit is σL/U = 1.74 ×10¹¹ Pa s (= 0.055 bar a), and the corresponding value of coefficient B in Glen’s flow law (3) is: B = U/(La³) = 7.92 × 10⁻²⁴ Pa_-3 s^-1 (= 0.25 bar_-3 a^-1). This value will be discussed later.

4.2.2. Shape of the cavity

Diagrams of the cavity profile have been made using the photogaphs taken at time t₀ then after one and four revolutions of the machine (times t₁ and t₄ in Figure 4). Figure 11 shows the variations of ρ (measured to the millimetre) versus θ for the points on the cavity roof. A slight increase of the cavity volume can be noted between the extreme measurements. The computed cavity corresponds well to observations, especially for the one obtained after one revolution (t = t₁ ).

Fig. 11. Radius ρ(θ) of the cavity roof: (A) cavity observed at t = t₀ (beginning of observations), (B) cavity observed at t = t₁ (after one revolution), (C) cavity observed at t = t₄ (after four revolutions). Solid line: cavity computed with n = 3; dotted line: polar radius of the bump profile.

4.2.3. Deformation of the marker network

The trajectories are solutions of the differential system:

(6)

Each trajectory is defined by its starting point M ₀ on the radius θ = π. The solution to Equations (6) is found starting from M ₀ at time t = 0 = t₀ and using the Runge-Kutta method. At time t, the current point M belongs to element e. The position of M at t + Δt is calculated using the nodal velocities solution to finite-element system (5) and the interpolation functions specific to element e. The time step Δt is taken small enough so that, if M crossses the boundary of e during Δt, its new position can be kept. This computation was performed for 19 stream lines, equally spaced at 5 mm on the radius θ = п. Figure 12 shows the deformation of the material line θ = 0, initially radial, after one and four revolutions.

Fig. 12. Computed deformation (n = 3) of the radial line θ = 0: (0) initial position; (1) after one revolution; (4) after four revolutions (this line is not part of the network which is already deformed at time t₀).

During the experiment, the ice trajectories are observed by following the displacement of each marker.

The comparison of the calculated and observed deformations of the network was made after one and four revolutions. The initial position of the network (instant t₀ in Figure 4) is shown in Figure 13. In order to reduce the computational burden, each individual marker trajectory was linearly interpolated between the two nearest among the 19 already computed trajectories. The results are shown in Figures 14 and 15. After one revolution the computed and observed networks are practically the same. After four revolutions the results remain quite good: the mean deviation in the angular position θ of the markers is 1°.

Fig. 13. initial position of the marker network (time t₀ in Figure 4).

4.2.4. Time description of the trajectories

More detailed checking of the flow can be carried out by studying how a marker moves along its trajectory with time. The study is restricted to the description of half a trajectory (0 > θ > π). About 50 photographs are available to measure the angular position θ(t) of a marker. The computed values are 12 times more numerous. The results for a marker on the innermost stream line are shown in Figure 16. As the differences between the observed and calculated θ are very slight after half a revolution, instead of θ(t) this figure shows the angular difference θ(t) to a point which would turn at the constant speed of rotation Ω of the machine, i.e. θ(t) = Ωt − θ(t). The comparison is made by taking the rejoining point on a bump as the starting point of the trajectory. It can be noted that when passing over the cavity, the differences θ are slightly greater for the marker than for the theoretical material point, which indicates a smaller increase in the observed

Fig. 16. Angular difference θ(t) = Ωt – θ(t) for a marker moving along the ice-bump contact and the cavity roof, o: observed; solid line: computed (n = 3) values.

velocity than that calculated with n = 3 in Glen’s flow law.

5. Discussion

5.1. Experimental conditions

Before discussing the results, it is necessary to look at how the assumptions used for the computation were handled during the experiment.

5.1.1. Steady-state conditions

The cold-room temperature is very difficult to maintain at a constant value over long periods. These temperature variations give rise to heat exchanges with the ice and are largely responsible for the torque variations (see Fig. 4).

On the other hand, the steady-flow assumption supposes that the shape of the cavities does not change. The measured increase in volume of a cavity between zero and four revolutions was about 100cm³, corresponding to a 1% decrease in the volume of ice, half of this variation taking place during the first revolution (volume of the cavity after four revolutions: 680 cm³). This decrease in the volume of ice was found to be mainly caused by melting at the bump contact, as heat was provided by the outer environment, further densification of the ice after time t = t ₀ being undetectable (Meyssonnier, 1983). Two relatively important melting phases correspond to the decreases in torque recorded at times t = t ₀ and t = t ₂ in Figure 4.

Finally, it must be noted that the polar angle of the rejoining point varied little from the beginning to the end of the observations, passing from 44° to 43.5°.

5.1.2. Boundary conditions

The boundary conditions used for the computation and for which experimental validity must be verified are:

• (a) frictionless ice-bump contact: in order to ensure the existence of a lubricating water film at the ice-bump interface, the melting temperature T _m must be below the ambient temperature. According to Lliboutry (1976): T_m = T₀ − βp_m − AN/w, where T₀ is the melting temperature at normal atmospheric pressure, p_m is the pressure at the ice-

  water interface, β ≈ 0.1°C MPa⁻¹ (Harrison, 1972), N is the number of saline equivalents per kg of ice, w is the water content, and A = 1.86°C kg mol⁻¹. With p _m = 6.5 Mpa, the mean normal stress corresponding to the measured torque, the drop in T_m is 0.65 C To reach the cold-room temperature N/w should be near 0.2 mol kg⁻¹. This is quite a probable value as w was probably very low. During the experiment no disappearance of air bubbles was noticed which would have indicated melting at grain boundaries (Barnes and others, 1971) and which was observed by Brepson (1979), who obtained a layer of blue ice totally free of bubbles when working at 0°C. Apart from the two short periods when melting was ascertained, ensuring perfect sliding, the friction should have remained very low at a temperature close to −1°C.

• (b) driving condition: the adherence condition on the driving ring seems to be well represented if one refers to the deformation of the marker network after one and four revolutions (see Figs 14 and 15). A computation from a separate finite-element model (Meyssonnier, 1983) indicated that the sliding velocity over a driving cog, corresponding to a mean drag of 0.07 Mpa, is negligible (u _b > 0.1 cm a⁻¹). However, it can reach 2 m a⁻¹ in zones where the tangential stress reaches its maximum value of 1 MPa (see Fig. 8).

• (c) plane-flow hypothesis: the observations made on pieces of the ice sample collected at the end of the experiment did not show any appreciable inclination of the markers from the vertical.

Fig. 14. Deformation of the marker network after one revolution (time t₁ in Figure 4): о observed: Δ computed (n = 3).

Fig. 15. Deformation of the marker network after four revolutions (time t₄ in Figure 4): о observed; Δ computed (n = 3).

5.2. Aptitude of the model for describing the observed behaviour

5.2.1. Value of exponent n in Glen’s flow law

The numerical results shown above were obtained with n= 3. This value is generally accepted for effective shear stress τ between 0.05 and 0.5 MPa (Barnes and others, 1971; Duval, 1976).

For higher values of τ, higher values of n have been reported in the literature (Steinemann, 1958; Dillon and Andersland, 1967; Barnes and others, 1971). In general, values close to n = 5 have been quoted. Goodman and others (1981) explained these high values by the fact that at high stresses (τ > 1 MPa) the dislocation sliding velocity is no longer a linear function of stress. According to Shôji and Higashi (1978) or Jones (1982), creep at high stresses is accompanied by generalized micro-cracking, initiated by the piling up of dislocations at grain boundaries.

In the present model, τ = 1 MPa with n = 3 corresponds to a dimensionless value τ = 1.18(U/BL)^1/3 or = 1.63(U/L). The extent of the region > 1.5(U/L) covers all of the part situated up-stream of the bump (see Fig. 8). In order to take into account the results of the literature concerning high stresses, the computations were run again with values of n higher than 3.

First, n was taken equal to 4. The corresponding map of does not differ considerably from that obtained with n = 3. Therefore, the comparison is made primarily on the total deformations. Figure 17, which represents the computed deformation of the radial line θ = π/2 after four revolutions, shows that this deformation is of course more important when n = 4, but that the differences are only appreciable near the innermost stream line. For a point on this stream line, the angular difference between its positions computed with n = 3 and n = 4 after four revolutions is 5.6°. As for the deformation of the marker array, the difference in behaviour, already noticeable after one revolution, is very clear in the final position. The model with n = 3 already has a slight tendency to overestimate the deformations and this tendency becomes very pronounced with n = 4. This value must therefore be discarded.

Fig. 17. Deformation of the radial line θ = π/2 after four revolutions, computed with n = 3 and n = 4.

In order to complete this work, the influence of a greater increase in n, but limited to the high effective shear-stress zones, was studied. To take such behaviour into account, Lliboutry (1969) and Colbeck and Evans (1973) have suggested adopting the polynomial relation = B ₁ τ + B ₃ τ ³ + B ₅ τ ⁵, in which the linear term represents the behaviour at low stresses. In this work we preferred to use Glen’s flow law in which n is a (step) function of т. The viscosity is still given as n = B ⁻¹ τ ¹⁻ⁿ where n = 3 and B = B ₃ for τ > τ ₅; n = 5 and B = B ₅ for τ ₅ ≤ τ; the coefficients B _j being taken to preserve the continuity of the viscosity. In practice, the exact non-dimensional value of the threshold τ ₅ which corresponds to a given stress cannot be fixed a priori (because the stress scale can only be fixed after calculation of the torque derived from the finite-element computation, and comparison to the measured one). Therefore, an approximate dimensionless value of τ ₅ is chosen, based on the results obtained with n = 3, and it is verified a posteriori that the actual value of τ ₅ corresponding to the experimental conditions is not too far from the desired threshold.

The computation was made for a threshold τ ₅ close to 1 MPa. The n = 5 region is then approximately limited by the contour = 1.5(U/L) in Figure 8. The ratio of the observed to the computed torques gives a pressure unit σ = 0.89 MPa and the corresponding actual threshold is τ ₅ = 1.05 MPa. For a marker in contact with the bump, the angular difference with its position computed with и = 3 after four revolutions lies between 2.5° and 3°, depending on its initial position. The difference between the final deformations of the marker network is less pronounced than in the case n = 4.

5.2.2. Value of coefficient B in Glen’s flow law

The torque-measuring system is very accurate but part of the measured torque is due to the friction of the mobile upper plate on its inner fixed support. This friction torque depends very much on the tightening and lubrication of the sliding pad. It can only be estimated by destroying the ice sample, allowing the ice to melt on contact with the bumps (this is done by raising the cold-room temperature and letting the machine turn without dismantling the upper plate). This procedure was not adopted here, in order to preserve the sample. A friction torque of 600 N m, known to within 400 N m, was assumed. This corresponds to the average of values obtained in other tests at the end of which the ice was melted. The resulting uncertainty in the measurement of torque is 10%. It remains acceptable for the estimation of the stresses involved, but the accuracy of S, computed for n = 3, is only of the order of 30%.

In other experiments at the end of which friction was measured, B was found to lie between 6.34 × 10_-24 and 1.27 × 10_-23 Pa⁻³S^-1 (0.2−0.4 bar_-3 a^-1).

The values of B computed with n different from 3 are:

(a) with n = 4: B = 7.6 × 10⁻³⁰Pa⁻⁴s⁻¹ (0.024 bar⁻⁴a⁻¹)

(b) with n = 3 if τ > τ ₅, n = 5 if τ ≥ τ ₅ = 1.05 Mpa;

B ₃ = 7 × 10⁻²⁴Pa⁻³s⁻¹ (= 0.22 bar⁻³a⁻¹

B ₅ = 6.3 × 10⁻³⁶Pa⁻⁵s⁻¹ (= 0.002 bar⁻⁵a⁻¹

The values corresponding to n = 3 are higher than Duval’s (1976), i.e. B = 3.17 × 10⁻²⁴Pa⁻³s⁻¹ (= 0.1 bar⁻³a⁻¹) at −0.05°C, for a low water content (w = 0.03%), and in a secondary creep state. This is due to the fact that at present deformation level tertiary creep is reached: re-crystallization and fabric formation contribute to increase the value of B. Comparison with Duval’s (1981) value for tertiary creep and zero water content, B = 6.3 × 10⁻²⁴Pa⁻³ s⁻¹ (= 0.2 bar⁻³a⁻¹ indicates that w remained very low during our experiment.

5.2.3. Steady-creep hypothesis

During the experiment the ice is obviously not in a state of steady creep: stresses at a given material point are very different up-stream and down-stream of the bump (see Fig. 8).

Figures 18 and 19 show the mean angular differences θ (as defined in section 4.2.4) observed for the markers which travelled close to the ice-bump interface. These averages are calculated as follows:

• (a) for each marker i, the observation allows θ at time t_n = nΔt to be calculated where Δt is the interval between. two photographs, Ω is the machine angular velocity, and θ ⁱ (t_n ) is the polar angle of marker ; measured at time t_n;

• (b) a value of θ ⁱ (θ_m) for is linearly interpolated through

• (c) assuming that θ ⁱ (θ₀) at θ₀ = π (see Fig. 18) or θ₀ = π/4 (see Fig. 19) for all the observed markers i, it is then possible to calculate an average ϴ(θ) corresponding to an Eulerian description of the flow.

Fig. 18. Mean angular difference ϴ(θ) for the markers along the cavity roof, compared to those computed with n = 1, 3, or 4 (θ = 0 is assumed at the bump summit for each marker).

Fig. 19. Mean angular difference ϴ(θ) for the markers along the bump contact, compared to those computed with n = 1, 3, or 4 (θ = 0 is assumed at point θ = π/4 for each marker).

Figure 18 shows ϴ(θ) observed for markers which have passed the bump summit and are moving along the cavity roof. A comparison is made with the values calculated for ice obeying Glen’s flow law with n = 1, 3, or 4. Figure 19 gives ϴ(θ) for markers which have passed through the rejoining point and are travelling along the ice-bump interface.

First of all, it must be noted that a purely Newtonian model is not suitable for simulating the ice behaviour. Along the first three-quarters of the cavity roof, the points calculated with n = 3 or n = 4 are ahead of the observed average (see Fig. 18). This tendency is reversed when approaching the next bump. For a marker in contact with the bump, the accordance with the computation for n = 3 is excellent, except on approaching the bump summit (see Fig. 19).

The difference in behaviour observed between a theoretical Glen body and actual ice, down-stream of the bump, could be explained by the influence of transient creep occurring on unloading (Duval, 1976; Duval and others, 1983). As for (re-)loading, no appreciable influence of Andrade creep (which appears upon loading of a virgin or annealed sample) can be noted. This result is in accordance with Lliboutry’s (1975) estimations.

6. Conclusion

The analysis of the observed and computed final deformations of the marker array allows for the assumption of a Glen’s flow law exponent n higher than 3 to be discarded, and all the more so since the actual deformation is higher than or equal to that which would take place if the ice-bump contact were perfectly frictionless.

However, it remains possible that n could take a value higher than 3 for effective shear stresses τ greater than 1 MPa: the region involved would be too small for a

significant difference in behaviour to be detected by our method.

A more detailed study on how the observation markers move along their trajectories shows the difference in behaviour between ice and a purely viscous body, and the influence of transient creep. It appears that only transient creep on unloading may be significant.

The comparison betweeen the angular velocities θ deduced from the delay curves in Figures 18 and 19, and the computed values (shown in Figure 20), indicates that these transient effects can be ignored when the basal flow velocity of a glacier in steady-state conditions is to be simulated.

Fig. 20. Deviation from the constant speed of rotation Ω of the angular velocity of markers on the inner stream line: values of Comparison between computed values and those derived from the observed values of ϴ(θ)

Acknowledgements

This study was carried out at the Laboratoire de Glaciologie et Géophysique de l’Environnement (C.N.R.S.-Grenoble, France), under the supervision of Professor L. Lliboutry, and with the financial aid of the Centre National de la Recherche Scientifique. Thanks are extended to Dr R. Brepson and Dr P. Duval, and to the late Professor P. Le Roy, for their helpful discussions and encouragement during the course of this work.