Preliminary investigations of centrifuge modelling of polycrystalline-ice deformation

Abstract

Glen’s flow law is a well-established general law for steady-state glacier ice deformation, and many laboratory tests and field measurements have been undertaken which have shown the generality of the law to be correct. In Nature, ice deformation is the response of the glacier/ice sheet to the applied self-weight stress of the ice mass (i.e. ice thickness, gravity and ice density) which produces a stress gradient within the ice column. Detailed experimental analyses of ice samples in the laboratory have until now only been undertaken using uniform stress fields in uniaxial or triaxial tests. Obviously the best method for investigating ice in the laboratory would be if stress gradients similar to those found in Nature could be replicated. In the following paper we describe the physical modelling of two (laboratory-prepared) isotropic, polycrystalline ice models (0.75 × 0.25 × 0.18 m) at enhanced gravity levels (80g) in a geotechnical beam centrifuge. Steel plate was placed on top of the ice model, replicating an overburden of approximately 36 m of ice (at 80g). Thus we were able to model the deformation of the lower 14 m of an ice mass approximately 50 m thick. Models are confined laterally by the Perspex strongbox walls, preventing lateral extension within the sample during testing. Models are unconfined on their downslope ends, rendering longitudinal stresses negligible. Deformation can therefore be treated as simple shear. Samples are instrumented with displacement markers and thermocouples. Values for A and n in the flow law derived from the experiments are reasonable and indicate the potential of this method for ice-deformation studies.

1. Introduction

A physical and theoretical understanding of the rheological behaviour of ice has long been the quest of glaciologists. It is fundamental to understanding how glaciers and ice sheets flow, and represents a major component of ice velocity in many ice masses. It is of great importance that we can model accurately the flow of glaciers and ice sheets and smaller ice bodies of the world, through the past, to the present and into the future.

At present, Glen’s flow law provides an adequate description of steady-state ice deformation for most applications (Glen, 1955). It relates the strain rate, thus:

(1)

where A is a constant related to temperature, grain morphology, liquid H₂O content and hydrostatic head, and has values between 10^–14 and 10^–21 s^–1 kPa^–3 (Paterson, 1994), τ_xz is the longditudinal component of the shear stress (where z is perpendicular and × is parallel to the surface slope) and n is a constant whose value is approximately 3.0 with a range of approximately 2–4 (Weertman, 1973).

However, transient creep is a significant factor both in laboratory tests which are designed to parameterize Equation (1) and in small ice bodies such as cirque glaciers and active rock glaciers whose ice residence time is of the same order of magnitude as the transient effects. Thus, a more general description of the rheology would combine Equation (1) with a transient creep component, ε_T, thus:

(2)

where ε_E is an instantaneous elastic response and is the steady-state (secondary) creep rate.

The transient component can be modelled using a term in which decays to zero at t = ∞ (viscoelastic) of the form:

(3)

where B = (1 + v)/2E _Young and τ_R is a relaxation time which depends upon various material properties. Another way of describing the transient creep component uses a power-law dependency (Azizi and Whalley, 1995) of the form

(4)

where A and m are creep parameters, A being the same as that used in Equation (1).

Clearly, where m = 0, steady-state conditions apply, and where m = 1, accelerating (tertiary) creep is assumed.

Such transient effects dominate the behaviour of small ice bodies in their initial stages. Reconstructions of the dynamics of small ice bodies such as alpine rock glaciers (Kääb, 1996) have used steady-state rheologies for simplicity. However, the importance of the transient creep effects has been demonstrated by Azizi and Whalley (1995), although parameterization of such effects remains problematic. Field sociologists have known for a long time that ice character varies markedly in real ice masses, especially within the lower ˚basal ice" layers of glaciers (Kamb and LaChapelle, 1964; Lawson, 1979; Knight, 1997) and within the deforming zones of alpine rock glaciers (Vonder Mühll and Holub, 1992). Intuitively, rheological variability must occur as a result of the differences in the physical character of ice found throughout ice masses but most profoundly within the lower basal ice layers. Experimental and field observations have supported this (Swinzow, 1962; Kamb and LaChapelle, 1964; Kamb and Shreve, 1966; Goughnour and Andersland, 1968; Hooke and others, 1972; Echelmeyer and Wang, 1987; Paterson, 1994). These rheological variations have not yet been incorporated into any glacier flow models, and there are at least three major reasons for this:

Relationships between ice character and rheology are difficult to investigate. Controls over rheology in ice masses are multivariate, and isolating the roles of individual factors has proved difficult in both field- and laboratory-based studies.

Incorporation of variability in ice character, and therefore rheology, into models of ice-mass motion has proven difficult due to lack of evidence relating to the systematic spatial variability of these properties. The most obvious case for such a development involves the inclusion of a separate rheology for the debris-rich basal ice layer that lies at the base of most ice masses.

Parameters defining a single constitutive relation are often ˚tuned" on an ad hoc basis in individual studies, with the result that an acceptable description of field observations may be obtained. Moreover, such parameterization from the study of samples is inherently discretized, and the apparently continuous functions used to describe ice deformation are actually a piece-wise superposition of independently derived data.

Further work is required to improve the applicability of the flow law through better parameterization. At present, studies of ice deformation are conducted by the instrumentation of either a natural feature in situ or small samples considered representative of that feature exposed to an artificial stress regime (Barnes and others, 1971; Hooke and others, 1972; Baker, 1978; Duval and others, 1983; Jacka and Budd, 1991). The difficulty associated with field-based investigations (including the often harsh conditions under which such measurements are made) obviously points towards a laboratory-based approach being very useful and one which must be utilized. Until now, laboratory-based physical modelling has been limited to the application of uniform stress fields to samples, with stress being varied in incremental steps. This forces the extrapolation from discretized, small-scale, test samples through a piece-wise superposition, up to realistic scales. However, the stress gradient acting through an ice mass (referred to in engineering terms as prototype) is a function of gravitational potential as described by Newton’s second law (F = ma). The steady-state deformation response of ice is non-linear, following a power-law flow relation, whereby the rate of deformation is proportional to a value between the second and fourth powers of the applied stress (Equation (1)). The most effective way to physically model the non-linear aspects of ice deformation would be to apply a non-uniform stress field to a sample, and the continuous function described by Equation (2) can then be better parameterized to describe the deformation. In the prototype the stress gradient is the product of ice density, ice thickness and gravity and so can be replicated in a model ice sample by exposing the sample to an enhanced gravity field. This can be achieved by flying samples of ice in a geotechnical beam centrifuge. In the following sections we shall introduce the concept of centrifuge modelling, the scaling laws associated with it and the results from two tests on isotropic polycrystalline ice.

2. Geotechnical Centrifuges: Theory

The geotechnical beam centrifuge is a unique physical modelling tool, the "power"of which has been known to geotechnical engineers for some years, but is only now gaining a profile in Earth Sciences. In many earth surface phenomena the role of self-weight body forces is crucial in explaining observed behaviour, i.e. where in situ stresses change with depth and the response of the material may alter in response to these stress differences. As suggested above, the best way to model physically ice deformation is to replicate the self-weight stress gradient experienced in Nature within the model as a continuous function. A scale model fixed onto the end of a centrifuge arm which is then rotated ("flown") rapidly experiences ˚an inertial radial acceleration field which, as far as the model is concerned, feels like a gravitational acceleration field but many times stronger than Earth’s gravity field" (Taylor, 1995). The model experiences zero stress at the surface (in our models this is the surface of the overburden steel plate), but this increases with depth, producing a stress gradient that is directly comparable to that experienced in field prototypes.

2.1. Scaling

Certain scaling issues must be addressed for the various time-dependent and time-independent processes, in order to relate the model to the prototype.

2.11 Linear dimensions

If a scale model of the field prototype (Fig. 1) is subjected to an inertial acceleration field which is N times Earth’s gravity, then at model depth h_m, the vertical stress in the model, s, is:

(5)

The same stress occurs in the prototype where σ_zp = σ_zm, so the depth of the homologous point can be calculated as:

(6)

Hence, the scaling relationship for macroscopic linear dimensions is 1: N. Problems can potentially arise with understress being experienced towards the top surface of the model, and overstress at its base, due to the proportional relationship between radial acceleration (a) and distance from the axis of rotation (r) a = ω ² r, where ω is the angular rotational speed of the centrifuge. This error was calculated to be ≤ 2% for the model geometry used herein (Irving, 2000).

Fig. 1. The correspondence of inertial stresses generated in a centrifuge model to those in prototype (Nature).

2.1.2. Grain-size and temperature

Equation (6) is applicable over macroscopic length scales only. However, a continuous flow-law function treats the sample as both a continuum and (via A), a collection of par-tides. This is problematic in that the interactions between the ice grains are many and complex, occurring through processes which are both diffusive and translative. The latter processes will scale as 1:1. There is much uncertainty over the diffusive processes, but recent research estimates the relationship to be in the range 1: N² to 1: N³ (Irving, 2000). In this experiment:

D is a significant fraction of sample thickness (8 mm ≤ D ≤ 10 mm in a sample 180 mm thick), in order to satisfy the assumption of purely macroscopic (dislocation) interaction between grains

the temperature is below –3.5°C. An Arrhenius-type viscosity relationship is not assumed (Mellor and Testa, 1969).

However, empirically derived values for A are well constrained for power-law creep at this temperature (Irving, 1999). The tests were of a short enough duration, with small strains, that physico-chemical interactions leading to re-crystalhzation were negligible.

2.1.3. Creep

The discipline of centrifuge modelling is still relatively new, having only really developed in earnest since the early 1970s. It has been very much the realm of the geotechnical engineer, so the scaling laws have been developed only in areas which are of interest to engineers. Gold-regions centrifuge studies have been restricted to geotechnical tests on permafrost (Ketcham and Black, 1995; Harris and others, 2000) mainly related to thawing and consolidation processes. A limited amount of work has been done on the modelling of sea-ice-structure interactions, with the ice actually being formed during flight at enhanced g levels (e.g. Lovell and Scholfield, 1986). Centrifuge modelling of the creep of frozen soils has received little attention, and no data exist on the potential time-scaling of ice creep. Time-scaling for creep will lie somewhere between 1:1 and 1: N ³, and Smith (1995) suggests that, because stresses and temperatures are identical between the model and the prototype, scaling should be 1:1. This, however, ˚has not been experimentally demonstrated" (Smith, 1995). Time may scale by N ² if crystal sizes are reduced by N for a pure ice sample. However, a reduction in grain-size by a factor of N will lead to an alteration in the dominant creep mechanism from power-law to diffusion. This could then lead to the modelling of non-power-law creep. The data obtained from these tests should allow the first comprehensive conclusions of time-scaling for ˚pure" ice creep in centrifuge modelling. The time-scaling for creep is explored further below.

3. Method

The aim of this research project was to establish the potential for modelling polycrystalline-ice deformation using the geotechnical centrifuge. In order to minimize our centrifuge testing time (by maximizing strain rate) and to model ˚realistic" ice thicknesses, we decided to attempt to replicate the deformation in a 50 m thick glacier. To reduce the errors (from over- and understressing due to differences in r), models have been made thin and loaded with plate-steel overburden which, having a density 7.85 times that of ice, reduces the overall model thickness substantially. Thus, only the lower part of the glacier is modelled, which is acceptable since shear stress and hence strain are concentrated here as described by Equation (1).

3.1. Test-bed configuration

3.11 Test box

The test-bed design is shown in Figure 2 and described in detail in Irving (2000). A slope of 15° was calculated as the optimum for maximizing the strain rate while utilizing the full dimensions of the box. (A basal shear stress of 106 kPa was developed at the bottom of the model.) The two ice models were placed side by side in the strongbox. It was assumed that with the expected strains for such a short test duration there would be negligible mechanical coupling between the two samples. The outward-facing sides of the samples were confined by the Perspex walls of the test box. Deformation can therefore be treated as simple shear. Models were unconfined on their downslope end.

Fig. 2. The strongbox design, with cooling system and models in place.

3.1.2. Cooling system

To maintain thermal equilibrium during testing, three MEECH-ARTEX vortex tubes were used to cool the strongbox. Each vortex tube is capable of generating an air stream at –20°C with a flow rate of 220 L mirT¹. One tube was directed below the base-plate of each model, and the third was passed into the upper void space above the models. Models were insulated along the base and ends by 25 mm plywood, and the top enclosed by an insulating lid (Fig. 2).

However, heat removal by these air streams was not sufficient to balance the rate of heat transferral to the test-bed due to the motion of the centrifuge arm (175 km IT¹) through the air mass in the centrifuge pit. At best, temperature could be maintained at approximately –7°C for an initial period of approximately 30 min and then rose gradually during the remainder of the test to –1°C over a 5 hour flight. The test would then be paused for 12–18 hours to return the model to a lower temperature. Models remained on the arm in the centrifuge pit and were cooled by the vortex tubes.

3.1.3 Instrumentation

The instrumentation for each model comprised an array of K-type thermocouples, four columns of deformation markers in each sample and a video camera to monitor the time development of strain during the test. No data were collected from the video camera, due to technical difficulties, and thermal data were likewise not recorded for the final 5 hours of the test series.

3.2. Model construction

A metal base-plate was fixed into a wooden mould. Grains of polycrystalline ice with a size range of 8.0 mm < D < 10.0mm were sieved into the mould, with chilled water (at 0°C) added at regular intervals such that the surface was dry when the top of the mould was reached. The models were cooled in a domestic chest freezer using spot cooling from the base and ambient cooling via its exposed upper surface. In order to maximize data acquisition, two different models were constructed:

Model 1a

Construction was as described above, with the cooling process interrupted so that the talik, formed between the ice growing from the top and the ice growing from the bottom of the model, could be partially drained. This resulted in a three-unit stratigraphy. The lowest layer (0–5 cm) was bubble-poor ice, the middle layer (5–9 cm) was poorly-bonded ice with a high porosity, and the upper layer (9–18cm) was bubble-rich ice.

Model 1b

The cooling process was not interrupted and there was no visible evidence of large-scale heterogeneities within the model.

Both models were considered to be isotropic with regard to their crystal fabric. Once frozen, the vertical array of four K-type thermocouples (accuracy ±0.1 K) was drilled on the central line of each model (one-third distance downslope from its head) with the thermocouples mounted at 6 cm intervals. The four columns of 5 mm diameter plastic cylinders (each 5 mm in height) were inserted for displacement measurement. These were located on the outer face, next to the thermocouple array, and halfway and three-quarters of the way downslope.

3.3. Testing

Models were flown side by side in a radial acceleration field equivalent to 80g. The ice models were 0.18 m thick, with an overburden of 0.058 m of steel. Flying at 80g and following the scaling law given in Equation (6), this provides stress levels equivalent to those which would be experienced in a glacier approximately 50 m thick. The stress field in the ice replicated that in the lower 14 m of the glacier. Temperature data were collected at 30 s intervals and downloaded to a personal computer, but after 9.5 hours the logger failed, so continuous temperature data were recorded for only part of the test. Total testing time was 14.5 hours completed during four separate flights.

Models were then sectioned to expose the marker columns. A lower no-slip boundary was assumed and displacements were measured relative to this. Four of the columns were destroyed during this process. Samples were taken pre- and post-test for subsequent textural analysis.

4. Results

As continuous strain data were not recorded, we cannot quantify the rate of decay of transient creep. Maximum measured displacements at 12 cm from the bed are: for model la, 9.54 mm; for model lb along the centre line, 6 mm; and at the edge of the model in contact with the Perspex side-wall, 4.2 mm.

With reference to Figure 3, the final deformation is interpreted for each recovered column of markers.

Fig. 3. Deformation profiles from columns la-3 (a), lb-1 (b), lb-2 (c) and 3–3 (d) (solid curved line), along with the viscoelastic approximation (solid line) and the visco-plastic approximation (dashed line).

1a-3, centre line, halfway downslope

Increased overall displacement is attributed to the fact that during model construction the ice grains were allowed to warm and adopt a more spherical morphology. This would decrease the viscosity of the sample (A in Equation (1)). Furthermore, the drainage caused a higher bubble content in the upper two layers and, more specifically, a lower shear strength in the high-porosity middle unit.

1b-1, outer edge, halfway downslope

Reduced longditudinal displacement was recorded which is attributed to lateral strain. Whilst it was hoped that the Perspex test-box walls would apply a confining stress, in reality the cool, dry air circulating in the test box combined with the heat transferred through the walls to erode the outer face of the models, resulting in some lateral strain. Hence, a small but significant departure from the plane-strain approximation occurs at the outer model boundaries.

1b-2, centre line, one-third distance downslope

The position of this column may have been sufficiently near to the upslope face of the sample for strain to have had an upslope component. Again a departure from the plane-strain approximation has occurred.

1b-3, centre line, halfway downslope

These data are the best constrained, in that they are the least affected by any boundary conditions.

5. Analysis

The deformation data were converted to a shear strain making the following assumptions:

The system can be described using a plane-strain approximation. All strain is in the downslope direction

As the total strain is small and rotation of markers could not be measured, simple shear is assumed. Hence there is no normal component to the deformation.

Temperature, although cycled through a known range, is assumed to be constant. Whilst it is acknowledged that this is unsatisfactory, the range is such that pressure melting is not significant at its upper limit and significant internal stresses are not induced by thermal expansion.

The incremental shear strain, ε = tanβ (where /3 is the ice-surface slope), was calculated as the ratio of marker downslope displacement, d, to marker spacing, h:

(7)

This was converted to a strain rate, by dividing by the total time over which the models were subjected to an acceleration of 80g, t = 52 200 s. Some markers were disturbed during exhumation, so the displacement curves are not as smooth as would be expected. The expected strain for each displacement marker was then calculated using both Equations (2) and (4). No instantaneous strain component in Equation (2) was observed, and its effect is not considered in this analysis.

Equations (2) and (4) were parameterized using reasonable values for A, B, n, τ_R and m. A least-squares regression was performed to fit the generated curves from Equations (2) and (4) to the strain curves measured from the models. The predicted deformation from both of the models is plotted in Figure 3, with the viscoelastic transient component predicting greater deformation at depth than the visco-plastic transient component, while the reverse is the case in the upper regions of the models.

Whilst the fitted curves describe the observed deformation well in the lower regions of the models, the departure in the upper regions occurs because both equations are highly sensitive to shear stress. It is a goodness of fit in the lower regions which is most desirable. The best-fit parameters for each model for each descriptive equation are shown in Table 1.

Table 1. Best-fit parameters for A B, n, τ _R and mused in the viscoelastic and visco-plastic approximations (see Fig 3)

6. Discussion

The material properties of the the two models were similar except for the stratigraphic differences within model la. The effect of this layering was to reduce the relaxation time, τ_R, which would be expected given that a zone with a lower bulk density will store less elastic strain. However, the magnitude of τ_R in both models is several orders greater than other reported values (Jellinek and Brill, 1956) which are typically 100–1000 s for stresses and temperatures similar to those in the tests described above. It is suggested that this is due to scaling. If there is indeed a scaling relationship between the diffusive (and hence microscopic-scale) processes that lies at the upper end of the range suggested in section 2.1, then the relaxation time proposed in Table 1 for a model in the centrifuge is reasonable. Comparing the values obtained in Table 1 with the values obtained by Jelli-nek and Brill (1956) suggests an approximate scaling of 1: JV¹⁷ with time. Whilst the accuracy of this value is open to debate, it closes the range of probable values by a factor of N to somewhere between 1: N and 1: N³.

The issue of grain-size scaling is worth addressing briefly. If grain-size is reduced by a factor N, then Smith (1995) suggested that time would scale with 1: N². As mentioned previously Duval and others (1983) stated that a reduction in grain-size to 1 mm for polycrystalline-ice samples will result in diffusive processes becoming the dominant mechanism of creep. A reduction of crystal size by N in our samples could result in models creeping by diffusive rather than dislocation processes. Thus, for these initial tests crystal sizes were not scaled down. We have not incorporated a grain-size reduction into our calculations, but, as is implied by the treatment arising from Equations (3) and (4), there is a significant effect arising from the microscopic processes which would serve to skew the scaling from an expected 1:1 relationship for purely dislocational processes. This may well explain the intermediate value arrived at above.

7. Conclusions

The application of a stress gradient to a sample is intuitively the best method for modelling ice-deformation processes in the laboratory. From the evidence presented above it is clear that this can be successfully undertaken using ice models in a geotechnical beam centrifuge. Additional steel overburden can be added to the upper surface of models in order to increase the normal and basal shear stresses to a level that would be experienced at the bed of some valley glaciers or rock glaciers.

The results of this experiment are very encouraging, with strains being recorded by means of a series of displacement columns. Strains within the models were small, and tests would be improved if longer-duration flights could be undertaken. To this end, an improved cooling system would allow test durations well in excess of the relaxation times to be undertaken.

Scaling issues have been clarified somewhat, but can only be resolved with a better description of the microscopic-scale interactions within the models. Further testing with data describing the time development of the deformation will allow a more accurate parameterization of the transient-strain equations. This, in turn, will allow these effects to be subtracted from the total strain development, and a better parameterization of the steady-state creep with depth can then be undertaken.

Curve fitting of the total measured strain to total predicted strain using a viscoelastic and visco-plastic approximation provides values of A and n which are well within an order-of-magnitude accuracy (Table 1).

Increased instrumentation of models, including pressure cells, pore-water pressure transducers and continuous strain measurement, will better constrain the physical processes occurring.

With further investigation and better experimental technique, it is hoped that centrifuge modelling can contribute significantly to the study of glacier rheology