Application of a Flow Model to the Ice-divide Region of Devon Island Ice Cap, Canada

Abstract

The steady-state flow model of Reeh (1988) is applied to a flow line that starts at the highest point of the Devon Island ice cap, follows the surface crest for 7.6 km, and then runs down the slope for a further 3.7 km. The effects of bedrock undulations, divergence of the flow lines, the variation of temperature with depth, and a basal layer of “soft” ice-age ice are taken into account. A flow law with n = 3 and a value of A close to that of Paterson (1981) is used. Longitudinal stress variations are neglected so that shear stress is calculated by the usual formula. It is estimated that these calculated values may be in error by at most 30%. Depth profiles of effective shear stress, and of the components of velocity and normal strain-rate, are presented at selected points along the flow line. These illustrate the large variations that occur near an ice divide and over bedrock undulations of amplitude comparable with the mean ice thickness. The model gives good predictions of the surface profile and of longitudinal and transverse surface strain-rates measured at ten points along the flow line. Predicted depth profiles of horizontal and vertical velocity components are compared with those measured in a bore hole. Comparison is limited by the fact that the model works in ice equivalent, whereas about 20% of the ice column consists of firn with different rheological properties from ice. The vertical velocity prediction is good. However, the model does not reproduce well the shape of the horizontal velocity profile, although measured and calculated fluxes differ only slightly. Predicted annual-layer thicknesses are within 15% of the measured ones in the upper half of the ice column, which consists of ice deposited in the last 1000 years. Predicted thicknesses in older ice are too small and the discrepancy increases with depth. This might indicate increased precipitation or, more likely, a thinner ice cap in the climatic optimum. However, it could also result from the fact that the layer of “soft” ice has been thinning continuously since the end of the ice age, so that the ice cap has never been in a steady state.

Introduction

Paleo-climatic histories obtained from ice cores have stimulated interest in flow near ice divides. Several bore holes have been drilled in such regions; moreover, flow conditions near the divide affect the deepest ice, even if the bore hole is far down the slope.

The assumptions underlying the classical analyses of flow in a parallel-sided slab (Nye, 1957) and of the steady-state surface profile of an ice sheet on a horizontal bed (Vialov, 1958) break down near an ice divide. Raymond (1983) used finite-element modelling to study the general features of flow there for the case of isothermal ice on a horizontal bed deforming in plane strain. Paterson and Waddington (1984) adapted this model to take account of uneven bedrock, variation of temperature with depth, and deformation transverse to the direction of flow, in order to calculate the thicknesses of annual layers in Devon Island ice cap for comparison with thicknesses measured in a bore hole.

Bedrock undulations influence the distributions of stress and strain-rate and thus also of annual-layer thickness. Yosida (1964) and Budd (1970) used perturbation methods to study flow down a uniform slope with small harmonic undulations superimposed on it. Whillans and Johnsen (1983) eliminated a dubious assumption in Budd’s theory and tested the revised theory against data from the Byrd Station strain network in Antarctica. They assumed linear rheology. This eliminated the need for special consideration of the region near the ice divide, a need that arises because longitudinal stresses change the effective viscosity of the ice when the flow law is non-linear. These analyses were based on solutions of the biharmonic equation for a stress function. Reeh and others (1985) and Dahl-Jensen (1985) used perturbation methods and a non-linear flow law to study flow in the Dye 3 area, south Greenland. The model ice sheet was divided into layers to take account of changes in the effective viscosity due to temperature and other factors.

In this paper we use the model of Reeh (1988) to study the steady-state distributions of stress, strain-rate, and velocity, and the variation of annual-layer thickness with depth, in the divide region of an ice cap on an uneven bed. The flow line considered starts at the highest point of the ice cap, follows the surface crest for some distance, and then bends around to follow the direction of steepest slope at right-angles to the crest. Flow along a crest has not been studied before, although the flow line at Camp Century, north-west Greenland, and the initial part of the one at Dye 3, south Greenland, are of this type, as is the flow line through Crête, a suggested site for a new deep bore hole in central Greenland. We apply the model to the Devon Island ice cap, because an unusual amount of field data is available for checking the predictions. First-order perturbation models, such as those mentioned above, could not be applied here because the amplitude of the bedrock undulations is comparable with the mean ice thickness. Moreover, to study flow along a crest by the finite-element method would require a complex three-dimensional model. Paterson and Waddington (1984) had to assume that the flow line started on the crest directly up-slope from the bore hole. Their results are therefore not strictly comparable with ours.

Field Data

The ice cap lies on the eastern part of Devon Island between lat. 74.5° and 75.8°N., and long. 79.5° and 85°W. It has an area of about 15 600 km² and a maximum elevation of 1900 m. Measured ice thicknesses vary from 200 to 1000 m. The ice cap has an east—west summit ridge, which overlies a region of high bedrock, and is drained by valley glaciers that descend to sea-level on its northern, eastern, and southern sides.

The study area extends along the ridge about 10 km westward from the highest point and about 3.5 km down each slope. Figure I shows surface contours; elevations were determined by optical levelling. Relative elevations are accurate to ±50 mm, although absolute values may be in error by as much as 25 m. Ice thickness was measured by radio echo-sounding with a precision of about ±5 m. Figure 2 shows bedrock contours.

Fig.1. Contours of the ice-cap surface in metres above sea-level, at 10 m intervals, with locations of levelling lines (broken lines), surface-strain measurements (solid circles), and the bore hole. The flow line studied is shown by a solid line.

Fig.2. Bedrock contours in metres above sea-level at 50 m intervals, with radio echo-sounding lines (broken) and the flow line (solid).

The flow line chosen for modelling passed through a bore hole to bedrock (“Hole 72”, 299 m deep). The line was drawn perpendicular to the surface contours. It follows the ridge for 7.8 km and then curves around and runs down the northern slope ( Fig.1). It coincides with surface-levelling and radio echo-sounding lines along the ridge and is always within 500 m of the nearest line on the slope. In addition, spot soundings were made every 10 m along the flow line for 300 m up-stream from the bore hole. Figure 3a shows bedrock and ice-equivalent surface profiles along the flow line. Ice-equivalent thickness, which is 17 m less than the true value, was used throughout the calculations. The curvature of each surface contour at its intersection with the flow line (Fig. 3b), required for calculating the convergence or divergence of the flow, was measured on Figure I. At the bore hole, the temperature at 10 m depth is −23 C and the basal temperature is −18.5°C. The ice can therefore be assumed to be frozen to the bed throughout the study area.

Fig.3. a. Bedrock and ice-equivalent surface profiles along the flow line. Vertical lines mark places for which profiles of stress, normal strain-rale, and velocity are presented (Figs 9–12).

b. Curvature of surface contours along the flow line.

c. Thickness of the layer of “soft” ice-age ice. as a fraction of total thickness.

Other input data needed are accumulation rate, ice temperature, and flow-law parameters. Accumulation rate was taken as 0.26 m/a, ice equivalent, at the highest point, decreasing to 0.23 m/a where the flow line leaves the ridge, and remaining constant thereafter. These values are consistent with the measured values which are numerous within 1 km of the bore hole but sparse elsewhere. This left some scope for tuning the model by adjusting the distribution of accumulation rate so that the flux at the bore hole was close to the measured value. The temperature T at depth Z was calculated from the empirical formula

(1)

Here H is ice-equivalent thickness and G = 0.025K/m. The surface temperature T _s was taken as −23.2°C along the ridge, increasing by 1K per 100 m decrease of elevation down the slope. Equation (I), due to Raymond (1983), predicts temperatures measured in the bore hole to within 0.2K.

The flow-law parameters can be determined if we use both surface and bedrock elevations as input to the model. The quantity u _m/HC _x can then be calculated at intervals along the flow line. Here u _m is the horizontal velocity component averaged over the ice thickness H, and C _x is given by Equation (24) of Reeh (1988). As Equation (22) of that paper shows, a double-logarithmic plot of u _m/HC _x against basal shear stress has a slope of n and an intercept proportional to A _r. Here, n is the flow-law index and A _r. the flow-law multiplier at a reference temperature. The flow law is thus

(2)

where έ ij, σ’ ¡j are strain-rate and stress-deviator components, τ is the effective shear stress, and F(T) is a function of temperature. Figure 4 is such a plot for points at intervals of 100 m along the flow line. A least-squares fit gives n = 2.9 ± 0.1 and A _r = 18.9 × 10⁻¹⁸ Pa^−2.9 a⁻¹ for T = −23°C If one takes the commonly accepted value n = 3 and keeps the centre point of the distribution fixed, A _r becomes 1.6 × 10⁻¹⁶ kPa⁻³ s⁻¹, which is close to the value of 1.3 × 10⁻¹⁶ kPa⁻³ s⁻¹ recommended by Paterson (1981, p. 39). We used n = 3, A ₋₂₃ = 1.6 × 10⁻¹⁶ kPa⁻³ s⁻¹, and an activation energy of 60 kJ/mol. The left-hand point in Figure 4 refers to the point 100 m from the start of the flow line. The fact that it lies far above the regression line may be evidence for transient creep at the low total strain near the summit of the ice cap.

This method of determining flow-law parameters is an improvement over that used, for example, by Budd and Smith (1981). Those authors omitted the factor C _x , which takes into account the variations in the shape of the horizontal velocity profile along the flow line, and thereby assumed that that profile has the “laminar flow” shape everywhere. Omission of C _x makes a major difference in our case; the value of n is reduced to 1.9.

Measurements of closure and tilting rate of bore holes in Arctic Canada and Greenland show that ice deposited during the last glaciation deforms two to four times more readily than ice deposited since (Paterson, 1977; Dahl-Jensen, 1985; Fisher and Koerner, 1986). A layer of “soft” basal ice, with a value of A _r four times that of the ice above, was therefore incorporated in the model. The thickness of this layer is known only at the bore hole. At other points, the value was obtained by calculating the depth to the 10 700 a isochrone (the date of the end of the glaciation) and scaling the corresponding thickness by the ratio of observed to calculated thickness at the bore hole. This ratio was about one-third. Figure 3c shows the fractional thickness of the “soft” layer at each point.

Measurements available for checking the model predictions are: (1) surface strain-rate components at ten points along the flow line, and (2) depth profiles of horizontal and vertical velocity components at the bore hole. Paterson and Waddington (1984) published a horizontal velocity profile based on a preliminary analysis of the inclinometer data. We subsequently found an error in the depths of the final set of measurements so the profile presented here is slightly different. A detailed discussion of the inclinometer data is in preparation. The vertical velocity profile in this paper also differs slightly from that given by Paterson and Waddington (1984). Both were obtained by differentiating the vertical velocities measured in the bore hole (Paterson, 1976) to give vertical strain-rates, subtracting the strain-rate due to firn compaction, and then re-integrating these strain-rates to give an “ice-equivalent vertical velocity”. In the 1984 profile, the integration was started at the surface using the measured total shortening of the bore hole in 1 year, converted to ice equivalent, as the vertical velocity. This conversion is uncertain because it is difficult to make representative density measurements in the surface snow. In the profile published here, the integration was started from zero velocity at the base, where the temperature is −18.5°C. Moreover, in extrapolating the lowest measured strain-rate to the base we have now allowed for an enhancement factor of 4 in the lowest 5 m of ice; this was not done in the 1984 analysis.

Results and Discussion

Tests of model predictions

The reliability of all the predictions depends on the validity of two basic assumptions: (1) the flow line is perpendicular to the surface contours, and (2) the direction of flow is independent of depth. We have no evidence on the first point. However, a theoretical analysis suggests that, over rough bedrock topography, the flow divide may be displaced, at least slightly, from the surface crest (Raymond, 1983). As regards the second point, the direction of the velocity vector changes by 7° over the lowest 100 m of the bore hole. This is unimportant. The change may be greater,

however, at some points on the ridge.

With flow-law parameters obtained as described in the preceding section, the model was then used in the alternative mode, namely to predict the surface profile (Reeh, 1988). This procedure is, in effect, an adjustment of surface slopes and ice thicknesses so that all points lie on the regression line in Figure 4. Measured and predicted profiles are shown in Figure 5. Discrepancies are small except in the last 2 km of the flow line where predicted elevations are too low. The reason is unclear; perhaps the ice in this area has developed a fabric favourable for shear so that the steady-state flux can be maintained with a reduced surface slope. Note that the surface slope has maxima over the tops of the bedrock hills at x = 4 and 8.5 km, as predicted by simple perfect-plasticity theory.

Fig.4. Plot used to determine flow-law parameters. The quantities are defined in the text. The best-fit regression line has slope n = 2.9. The line drawn, which is for n = 3. is barely distinguishable from it at the scale of the diagram.

Fig.5. Predicted (curve) and measured (points) ice-equivalent surface elevations along the flow line. The arrow shows where the flow line leaves the ridge.

Figure 6 shows measured and calculated strain-rate components along and perpendicular to the flow line. As expected, the calculated values show more scatter than the observed ones, which are averages over distances of about 450 m. In particular, the averaging can explain the apparent overestimation of the longitudinal component at the bore hole (x = 8.8 km) and immediately up-stream. The model underestimates the transverse component for values of x between 1 and 5 km, but otherwise the agreement between prediction and observation is satisfactory. Transverse strain predominates along the ridge. The flow lines diverge everywhere, except at x = 9 km and in the last 2 km, where they are approximately parallel. The longitudinal strain-rate is compressing at some places on the ridge. These correspond to saddles in the bedrock ridge beneath the crest of the ice cap ( Fig.2) where an exceptionally high proportion of ice flows off to the sides.

Fig.6. Predicted (curves) and measured (points) surface strain-rate components along the flow line. The arrow shows where the flow line leaves the ridge.

In Figure 7, predicted velocity profiles are compared with those measured in the bore hole. The break in slope of the predicted profiles occurs at the top of the “soft” basal layer of ice-age ice. One feature must be emphasized. In the model, all the material is assumed to be ice. In reality, the uppermost 65 m of the ice cap is firn and the density does not reach 900 kg/m³ until a depth of about 100 m. Because firn has different rheological properties from ice, the model predictions are not expected to agree with observation in the upper layers. For example, the “softness” of low-density firn explains the curvature of the horizontal velocity profile near the surface.

The shape of the predicted horizontal velocity profile differs significantly from that observed. The predicted velocity gradient (shear strain-rate) is too large in the upper half of the ice column and too small over much of the lower half. This might result from failure of a basic assumption in the model, namely that variations in longitudinal stress can be neglected so that shear stress is given by the simple formula τ_xz , = ρgzα. In fact, the greatest perturbations in shear stress are expected near the top of a bedrock hill, where the bore hole is located, and in bedrock hollows. We estimated the magnitude of the correction to τ _xz in the following way. The values of the normal stresses σ _x and σ _z , calculated by the model, were used to determine ∂(σ _x − σ _z )/∂x. This quantity was then integrated with respect to z to obtain a first correction to τ _xz . This procedure follows from equations given by Reeh (1988). The correction term for the bore-hole site is positive near the bed and negative above a relative height of 0.15. This is the right direction to explain the difference between the velocity profiles, but the correction is not large enough. An additional factor is that the ice below about 180 m has a single-maximum fabric with near-vertical alignment of c-axes; this becomes progressively stronger as depth increases (personal communication from R.M. Koerner). This fabric would make the measured shear strain-rates greater than those predicted by the model in which the ice is assumed to be isotropic.

In spite of the differences in shape between observed and predicted velocity profiles, the calculated flux (139 m²/a) is close to the measured one (132 m²/a, ice equivalent). However, this agreement has to some extent been forced by adjusting the mass balance, as explained previously.

The measured and predicted vertical velocities (Fig. 7b) can be compared only below the depth (relative height 0.7) where the density of the material reaches 900 kg/m³. The “measured” profile was derived from the vertical velocities measured in the bore hole, after correction for firn compaction by the method explained in the preceding section. Because firn has different mechanical properties from ice, however, the velocity distribution calculated in this way is not the same as would be found in an ice cap of the same ice-equivalent thickness but consisting entirely of ice. For example, firn compaction accounted for the total deformation observed above a depth of 27.5 m (Paterson, 1976). If the uppermost four points in Figure 7b are omitted, all the measured velocities are within one standard error of the values predicted at the bore hole. However, all but one of the measured values are less than predicted. The predicted values for x = 8.7 km, 100 m up-stream from the hore hole, are also shown; these are all less than the measurements. The prediction for x = 8.75 km would therefore fit the observations almost perfectly. We believe that this displacement is an artifact. In the model, surface and bed slopes were determined by differentiating cubic spline fits to the measured surface and bedrock profiles. The slopes calculated for x = 8.75 are closer to those measured at the bore hole than are those calculated for the bore-hole site itself.

Fig.7. Velocity components measured in the bore hole (points) compared with predictions for them (solid curve) and also the predicted vertical component for a point 100 m up-stream (broken curve). Relative height is in ice equivalent. The firn—ice transition is at a relative height of 0.83 and density remains constant below 0.70.

Fig.8. Predicted values of effective shear stress at the surface and base. The arrow shows where the flow line leaves the ridge.

The predicted horizontal velocity profile for x = 8.7 km has the same shape as the profile for x = 8.8 km. So the discrepancy between the measured and predicted horizontal velocity profiles cannot be explained by the above-mentioned displacement.

Distribution of stresses

Figure 8 shows the computed distribution of effective shear stress along the flow line. The base curve is effectively a curve of basal shear stress τ _xz , because shear predominates near the bed, shear between adjacent flow lines τ _xy is neglected, and the transverse velocity component is assumed to be zero at all depths so that τ _yz is also zero. The base shear stress ranges from 0 to 50 kPa along the ridge and averages about 80 kPa down the slope. Along the ridge, relatively high values occur where the bed slopes downwards in the direction of flow, whereas values are low where the bed slopes upwards. This suggests that ice in the hollows is relatively inactive. Normal stress deviators predominate near the surface with the transverse component larger than the longitudinal one along the ridge and smaller on the slope. However, the surface value of the effective shear stress does not vary much along the flow line. Figure 9 shows how effective shear stress varies with depth at selected points along the flow line. On the slope, the effective shear stress increases with depth as previous analyses have shown (Nye, 1957). Along the ridge, in contrast, the effective shear stress decreases with increase of depth as Raymond (1983) found at a symmetric ice divide on a horizontal bed.

These distributions of stress, and also those of strain-rate and velocity, of course depend on the basic assumption that variations in longitudinal stress can be neglected. We have checked the validity of this by calculating the correction to the value of basal shear stress obtained from the simple formula. (See preceding section.) At worst, the correction is one-third of the basic term at a few places on the ridge where dH/dx is large. At these places, neglect of the correction term tends to enhance the pecularities of the various profiles.

Distribution of normal strain-rates

Figure 10 shows the depth variation of the normal strain-rate components at various points along the flow line. The transverse (y) component is always zero at the bed because the transverse velocity component is zero.

The profiles for x = 0, 3, and 4 km are similar: a large transverse component έ _yy at the surface, decreasing steadily to zero at the bed, and a small longitudinal component έ _xx which may be either extending or compressing and may increase or decrease with depth. At x = 5.1 km, έ _xx increases steadily with depth. Here, the bed slopes steeply downwards in the flow direction and the large extending strain-rate near the bed can be explained by a simple geometrical argument. If u is the x-component of velocity, B(x) is the z-coordinate of the bed, and suffix b denotes the value at the bed

(3)

With z measured positive upwards, ∂u/ ∂z is always positive and dΒ/dx is negative on a down-slope. In this case, (έ _xx )_b is positive. This is the contribution of bed topography to the longitudinal strain-rate; it must be added to the positive value required to maintain steady state in an accumulation area. Thus έ _xx near the bed can be much greater than its surface value at places where both the bed slope and the shear strain-rate in the basal ice are large. The profile at x = 8.4 km, where the bed slopes steeply upwards, shows the opposite effect. Here, έ _xx changes sign and becomes compressing near the bed because dB/dx is now positive. As a result, the vertical component έ _zz becomes extending near the bed and this is how the ice can move over the hill. Comparison of the profiles for x = 7.6, 8.4, and 8.8 km shows the large changes in longitudinal strain-rate distribution as the ice moves over a bedrock hill. At x = 11.3, where the bed is approximately horizontal, the profile of έ _xx resembles that given by Raymond’s (1983) formula for “laminar flow”. Because this formula applies to isothermal ice, however, it underestimates strain-rates in the lower half of the ice column.

Fig.9. Variation of effective shear stress with depth at selected points on the flow line. The break in slope of each curve occurs at the top of the “soft” basal layer.

Fig.10. Variation of longitudinal, transverse, and vertical components of strain-rate with depth at selected points on the flow line. The break in slope occurs at the top of the “soft” basal layer.

These results show the wide variation in strain-rate profiles that can occur near an ice divide, along a ridge crest and, more especially, where the bed topography is rough. In these cases, time-scale for ice cores calculated from simple flow models are unlikely to be of much value. For instance, if the ice travels over an upward-sloping bed for an appreciable part of its path, the annual-layer thickness near the bed may even increase with depth. The results also show that, where the bed is rough, the assumption, commonly made in reducing inclinometer data, that the normal strain-rates at depth do not exceed the values measured at the surface, may break down.

Fig.11. Depth variation of horizontal velocity component, normalized by dividing by its depth-averaged value, at selected points on the flow line. The break in slope occurs at the top of the “soft” basal layer.

Distribution of velocity components

Figures 11 and 12 show velocity profiles at points along the flow line. Variations in shape of the horizontal velocity profiles can be discussed conveniently in terms of the parameter ϕ(1), the ratio of the surface velocity to the mean velocity. For “laminar flow” of isothermal ice with n = 3, ϕ(1) = 1.25. Higher values of ϕ(1) indicate that shear strain is more uniformly distributed over the ice thickness than in laminar flow, whereas values of ϕ(1) less than 1.25 show that shear strain is more concentrated in the basal layers.

Near the highest point of the ice cap, the profile of horizontal velocity u has an inflection at approximately mid-depth and ϕ(1) has a value near 2. In Raymond’s (1983) approximate analytical solution for the divide, u varies linearly with depth. Raymond argued, however, that the profile must have an inflection near the bed, and his numerical results at a distance of half the ice thickness from the divide show this. The inflection in our curve for x → 0 is much more marked. Our profiles at x = 3 and 7.6 km also have inflections. The surface slope, basal shear stress, and horizontal velocity along the ridge all have minima at these points and the flow lines diverge strongly. These points are, in effect, subsidiary ice divides, so it is not surprising that the velocity profiles resemble that near x = 0. The vertical velocity profiles at all three points are similar. The value of ϕ(1) varies along the ridge; values are high where the normal stress deviators have high values, especially near the surface, whereas ϕ(1) has low values where the basal shear stress predominates. Down the slope (x > 7.8 km) ϕ(1) decreases with increasing distance, in agreement with Raymond’s result.

Horizontal velocity profiles with high values of ϕ(1) are associated with vertical velocity profiles that are concave downwards, that is, vertical motion tends to be concentrated in the upper layers. This occurs at places on the ridge where έ _yy has high values near the surface and, on the slope, where έ _xx is high. Both of these result in high near-surface values of έ _zz . Raymond (1983) found that, if the bed is horizontal, The vertical velocity profiles are always concave downwards, although the concavity becomes less marked as distance from the divide increases. In the present case, the profiles at x = 5.1 and 8.8 km are nearly linear and that at x = 9.6 km is concave upwards. These are at places where the bed slopes steeply downwards so that έ_zz is strongly compressing there. At x = 8.4 km, the vertical velocity is directed upwards in the lowest 35% of the ice thickness so that the ice can flow over a bedrock hill.

Raymond (1983) suggested two reasons for variations in the shapes of the velocity profiles: (1) a wide range of longitudinal strain-rates which results in variations in the depth distribution of the effective viscosity, and (2) a wide range of longitudinal stress gradients which changes the depth variation of shear stress τ _xz . He stated that the second was the more important, although it is not clear how to assess the relative importance of two quantities of different dimensions. Our model cannot shed much light on this question because we have excluded the second effect by our assumption that τ _xz is always given by the standard formula. However, as stated above, our tests show that inclusion of the longitudinal stress-gradient term would change the value of τ_xz by, at most, 30%.

The inclusion of a “soft” basal layer in the model increases the amount of adjustment of the ice column to the rough bedrock topography that can take place in the basal ice. We therefore expect that, without this basal layer, the variations in velocity and strain-rate profiles along the flow line would be greater than shown here.

Annual-layer thickness

Annual-layer thickness as a function of depth in the bore hole was calculated as the reciprocal of the derivative, with respect to depth, of the time-scale (the age—depth relation). Particle paths and the time-scale were obtained by numerical integration of the velocity field. It was found that the ice 14 m above the bed originated about 1 km upstream. Calculated layer thicknesses are compared with the measured ones (personal communication from D.A. Fisher and R.M. Koerner) in Table I. Down to a depth of 136 m, where the ice is about 900 years old, calculated values are within 10% of the measurements, except for one point. Below this, the measured thicknesses exceed the calculated ones and their ratio increases with depth. At 267 m, for example, the measured layer thickness is six times the calculated one. As a result, the calculated age is about 11 000 a as against the true value, obtained from measured layer thicknesses, of 4800 a. Paterson and Waddington (1984) obtained a similar result with their finite-element model; calculated thicknesses were within 10% of the measured values down to 170 m (age 1300 a) with the ratio of measured to calculated increasing steadily below this. The flow line used in the finite-element model ran directly up the slope to the crest of the ridge; flow along the crest was ignored. This difference may account for the minor difference between the results. Paterson and Waddington concluded that the ice cap in the bore-hole region had been in a steady state with the present accumulation rate for about the last 1300 years and, before that, accumulation rate was greater by an unknown amount.

Fig.12. Depth variation of vertical velocity component at selected points on the flow line. The break in slope occurs at the top of the “soft” basal layer.

We believe that even these vague conclusions are doubtful. Application of a simplified non-steady-state flow model (to be presented elsewhere) indicates that the measured layer thicknesses are compatible with a gradual thickening of the ice cap near the bore hole from about 100 m 5000—6000 years ago to the present thickness of about 300 m. This compares well with the 250 m of thickening needed to reconcile the Devon Island and Camp Century oxygen-isotope profiles over the past 5000 a, on the assumption that the climatic change over this period was the same at both stations (Paterson and others, 1977). These ice-thickness changes are surprisingly large. Two other factors may help to explain why the annual layers in the deep ice are much thicker than predicted. First, the ice below 180 m, which has a single-maximum fabric with near-vertical alignment of c-axes, will deform less readily under vertical compression than the isotropic ice assumed in the model. The second factor is time variation of the strain-rate/depth distribution. Since the end of the glaciation, the layer of “soft” ice-age ice has thinned continuously and been replaced by an increasing thickness of “harder” ice on top. In the past, when the “soft” layer was thicker than now, it would have contributed a larger proportion of the total deformation needed to move the ice down-slope and over bedrock hills. Thus, a steady-state model, with a “soft” layer of its present thickness, will overestimate the total deformation of the “hard” ice.

Table.1. ANNUAL-LAYER THICKNESSES

The existence of two layers of different effective viscosity, whose relative thickness varies with time, means that the ice cap has never been in a steady state since the end of the glaciation (Reeh, 1985). To compensate for the steady increase in mean viscosity as the “soft” layer thinned the thickness and/or surface slope must have increased steadily to carry off the same accumulation rate. In addition, the increase in ice temperature associated with the end of the glaciation must have influenced the flow for some time. It follows that a steady-state model will not predict the correct layer thicknesses even if accumulation rate remains constant. We therefore support Paterson and Waddington’s view that this is a poor method for determining past accumulation rates. In addition, our study shows the severe limitations of steady-state flow models, even detailed models supported by large amounts of field data, for calculating time-scales.

On the positive side, our model, with well-established values of the flow-law parameter, successfully predicts the surface profile of the ice cap and the widely varying values of surface strain-rates along the flow line. It also shows the characteristics of flow along a gently sloping ridge and the wide variations in depth distribution of stress, strain-rate, and velocity produced by flow over bedrock hills and valleys.

Acknowledgements

The field data were obtained while the second author was with the Polar Continental Shelf Project, Department of Energy, Mines and Resources, Canada. A substantial part of the flow model was developed while the first author was with the Geophysical Institute, University of Copenhagen, and was sponsored by the Danish Commission for Scientific Research in Greenland and the Danish Natural Science Research Council.