Introduction

Relations describing the mechanical response of ice to an applied stress are fundamental for understanding ice-sheet dynamics. The commonly used relation is the Vialov (1958) differential equation, which links the derivative of the ice-deformation velocity with respect to depth in the ice sheet to the effective stress τ at power n, where n is called the Glen parameter, and to the ice temperature via an Arrhenius-type relation depending on an activation energy Q (i.e. of the type B(T) = B ₀ exp(−Q/RT)). Although the relation is well-supported by experiments or in situ measurements, a very large set of empirical values is found in the literature: values of the Glen parameter vary from 1 to 4.5, with a mean value of 3 (Paterson, 1981). Homer and Glen (1978) reviewed 23 values of activation energy between 40 and 135 kJ mol⁻¹, and found a new value of 75 kJ mol⁻¹.

Concerning the Glen parameter, a quasi-Newtonian viscosity (n ≈ 1) is suggested for polar ice as long as grain growth occurs (Pimienta and Duval, 1987). This creep concerns the first hundreds of m at the top of the ice sheets. A power law creep with n = 3 is expected when continuous recrystallization occurs; it is named tertiary creep. Results concerning the intermediate ice layers are uncertain: rotation recrystallization involving sub-grain-rotation and grain-boundary migration was suggested by Pimienta and Duval (189) for these kinds of ice. Several interpretations of field data suggest a flow with n = 3 (Paterson, 1983; Reeh and others, 1985) whereas others have given n from 1 to 2 (Lliboutry and Duval, 1985; Pimienta and Duval, 1987). Moreover, the exact boundary between different types of creep and their relative importance is not very well known. Creep associated with rotation recrystallization probably occurs for 80% of the ice sheets but its exact importance is difficult to evaluate. In particular we do not know whether it is significant in the deformation processes which occur mostly at the bottom.

Activation energy seems to depend on temperature. Creep tests made by Duval and Le Gac (1982) gave a value of 78 + 4 kJ mol⁻¹ for temperatures between −10 and −25°C. A larger value is expected for warmer temperatures. Finally, the pre-exponential factor (B₀ ) depends on crystal structure, on concentrations of impurity and other factors (Budd and Jacka, 1989).

Laboratory measurements of these parameters can be done easily only for high stress, because of the time needed to achieve measurable deformation. They thus usually yield large values of n (Lliboutry and Duval, 1985), in contrast to borehole studies which indicate a Newtonian flow.

Empirical studies from estimations of flow velocity and stress are a means to get information over a wide range of stress, strain rate or temperature. The major difficulty lies probably in the stress computation which requires a precise ice-sheet surface topography. This was attempted by Young and others (1989) and Hamley and others (1985) using the topography of Zwally and others (1983) and of Drewry (1983). Neither of these studies took into account the temperature dependence; both found a high value of n for the whole data set.

This paper re-estimates parameters of the constitutive flow law from the same empirical approach, using a very precise topography deduced by inversion of satellite altimeter data (Remy and others, 1989). We focus on the area between 1500 and 3000 m, where the topography is very precise. In addition, temperature dependence is taken into account. The first section of the paper is devoted to formulating the constitutive law: the formulation must be complete enough to be realistic and simple enough to allow regression. Then the stress, U/Η (ratio between mean balance velocity and thickness) and temperature are computed along ice flowlines derived from surface topography, to determine flow parameters.

Finally, rheological parameters n and Q are derived by least-squares regression and by assuming that, locally, they do not vary along a flowline. This semi-empirical approach is applied to four sectors of Antarctica (Fig. 1).

Constitutive Equation

In this section we write the constitutive equation linking the velocity at a given position x, along a flowline, to other ice-sheet parameters. The whole presentation is based on steady-state hypothesis.

In cold ice sheets, temperature and velocity are closely related through the Vialov (1958) differential equation which links the derivative of the “deformation velocity” u(z) with respect to height above ice bottom z, to temperature Τ and basal shear stress τ:

where Η is the ice thickness and R is the gas constant. T_m is the melting temperature at the bottom of the ice sheet (T_m = 273 – H/1503, where T _m is expressed in K). The stress is taken as:

where ρ is the ice density, g is the acceleration of gravity and ∝ is the surface slope along a flowline direction. We neglect other stress components. First, the stress is calculated over a distance of the order of 50 km: the effect of longitudinal stress and stresses due to the irregular bedrock topography are then smoothed (Young and others, 1989). Secondly, we avoid domes and coastal regions where shear stress is probably not dominant. To integrate Equation (1), we choose the formulation of Lliboutry (1979). He noted that maximum deformation occurs in the first 100 m above the bottom where temperature varies linearly with depth, and that the integration depends only on the bottom temperature gradient. He deduced the velocity profile as:

where

and U is the averaged velocity over the ice column such that

with

and

T_b is the mean bottom temperature averaged over the first 5% of the bottom ice (which corresponds to about the first 100 m) and G ₀ is the bottom temperature gradient, taken as 0.022°C m⁻¹ corresponding to Φ = 50 mW m². Note that (p + 2) arises in Equation (5), which is linked not only to the flow parameters n and k but also to the ice thickness H. The pre-exponential factor B ₀ is independent of the ice-column characteristics.

Computational Scheme and “Flow Parameters”

Assuming that they are constant, n, k (or equivalently Q) and B ₀ can be estimated by regression from Equation (5) if U/H, the basal shear stress τ and the bottom temperature T_b are known along a set of geographical positions. An iterative scheme will be used. First, ice flowlines are deduced from the direction of maximum surface slope of the ice topography (Fig. 1 ). Between two flowlines, the steady-state continuity equation is applied step by step:

where x is the flow direction of the ice divide to the coast and dx is the step, chosen as 20 km, the overbar designing mean value between x and x + dx. l(x) is the distance between two flowlines. Ice thickness H(x) is derived from Drewry (1983), the surface ice-sheet topography from Remy and others (1989) and accumulation rate b(x) from Radok and others (1987). U(x) and U/H can then be computed step by step.

Fig. 1. Topographic map of the selected sector of Antarctica, deduced from Seasat data, as explained in Remy and others (1989). Bold bolines are each 100 m. The figure shows also the four selected flowlines.

The flowlines are derived from the surface slope of the topography: the initial distance l(x) is 20 km. The integration starts from the divide line, using the topography of Drewry (1983). North of 72° S, our altimetric topography is used, and only the results corresponding to the latter are discussed. Surface slope ∝ is deduced from the topography by fitting a biquadratic form over 25 km around each point; stress is then deduced from Equation (2). The calculation is stopped at 68° S because the quality of the altimetric topography is seriously degraded north of this line.

The bottom temperature must be estimated from the thermodynamic equation (Lliboutry, 19877) which, in the steady-state hypothesis, reads

where k is the thermal diffusivity, C the specific heat capacity and w the vertical velocity. In this equation, we consider vertical diffusion (first term), horizontal advection (second term), vertical advection (third term) and dissipation (last term), w is calculated from mass continuity:

In Equations (9) and (10) u(z) is estimated as in Equation (3), from U computed at the same place by Equation (8). This has been proven a sufficient approximation to estimate bottom temperature (Ritz, 1987).

The boundary conditions are, at the base of the ice sheet, the bottom temperature gradient (G ₀), and at the surface the temperature T _s is derived from Radok and others (1987). ∂T/∂x at the surface is supposed to be equal to αλ where λ is the atmospheric vertical temperature gradient, assumed to be 0.0115 deg m⁻¹ in the selected region (Huybrechts and Oerlemans, 1988). The initial conditon is u(z) = 0 at the dome, corresponding to a purely vertical flow induced by accumulation at the surface.

Along each flowline, we then compute U/H, τ and T_b from Equations (5), (8), (2), (9) and (10). Because u(z) depends on p and k, T _b is calculated first by using initial values for these parameters (n = 3, k = 0.1; Radok and others, 1987). The whole process (i.e. calculation of U/H, τ and T _b up to re-estimation of n and k as discussed below) is then itereated with the new couple of values. It converges in two iterations and the results were shown to be unsensitive to the choice of the initial values. The calculations have been achieved along four different flowlines (Fig. 1). The values of the involved parameters along flowline Β are illustrated in Figure 2.

Rheological Parameters

In this section, we analyze the best-fit values for the activation energy Q (via k, Equation (6)), n, and the pre-exponential factor B₀ . We only consider data where T _b is less than −10°C. Thus we exclude areas where melting and strong basal sliding might occur. Figure 3 shows U/H versus stress for the four flowlines: the link between these parameters is not linear even in a logarithmic scale. Moreover, for flowlines A or D, one observes an increase of the velocity when stress remains constant. A line corresponding to n = 3 is shown for reference: as a matter of fact, if one does not take into account the temperature dependence (k = 0), one finds n between 2.7 and 3.3 for the four flowlines. This is in good agreement with the values of Young and others (1989) and Hamley and others (1985).

Fig. 2. Values of bedrock and surface topography (1), accumulation rate (2), surface temperature (3), stress (4), U/H (5) and mean bottom temperature over the bottom 5% of the ice layer (6), along flowline Β (Fig. 1).

Fig. 3. U\H values for four selected regions A, Β, C and D (Fig. 1) versus stress, on a logarithmic scale. If temperature dependence is not taken into account, n = 3 is found for the four regions. When the effect of temperature is considered, smaller values of n are found (see Table 1).

However, variations of stress and temperature change in sympathy along flowlines (Fig.2): a principal component analysis from the three-dimensional correlation matrix of these parameters (using log (U/H), log τ and T_b ) shows that the first eigenvalue represents more than 99.99% of the variances. As a consequence, the calculated values for k and n are strongly correlated.

The best-fit values for n and k for each flowline were derived by iterative least-squares calculation. This allows us to take into account the non-linear dependence of p on n and k (Equation (7)). As shown in Table 1, n is found between 0.64 and 1.27, and k between 0.10 and 0.12, which corresponds to an activation energy between 65 and 80 kJ mol⁻¹. This value is close to the value of 75 kJ mol⁻¹ of Homer and Glen (1978) and of 80 kJ mol⁻¹ of Duval and and Le Gac (1982).

The pre-exponential factor B ₀ shows wide variations depending on the selected flowline. Ice-sheet dynamical models also show that B₀ is very dependent on the flowline (Huybrechts and Oerlemans, 1988). However, comparison between our result and literature values is not direct because of the normalization by (p + 2) (Equation (5)).

Confidence Tests

We now discuss the different factors which may induce some errors in the derived parameters and the confidence in the above results. Because of the regular behaviour of U/H and stress, in the case of flowline B, the discussion is directed to this case first. Both ice path and accumulation rate have a cumulative effect along the flowlines, because of their contribution to the balance velocity. Their errors are propagated in all the downstream results. On the contrary, bottom temperature, stress and ice thickness only have a limited extent. Thus, a poor knowledge of the latter will introduce minimal error downstream.

In a first test, the stress is calculated over a scale of 100 km, rather than 50 km, only the n value is diminished (Table 1, case 1). A series of tests concern the computation of T_b (Table 1, cases 2–5). First, we replace the vertical veolocity profile (advection is the dominant vertical heat-transport process) derived from Equation (10) by a linear profile, set to the accumulation rate at the top of the ice sheet and to zero at the bottom. We also test the effect of the atmospheric vertical temperature gradient (λ) and values of geothermal flux Φ. In each case bottom temperature variations are within a few degrees and the derived rheological parameters are within 20% of the nominal ones. Case 2 suggests that the approximation of Equation (3) is acceptable.

The very dense network of airborne radio-echo sounding missions near region Β leads to a precision of the bedrock topography better than 50 m on the 100 km scale (Drewry, 1983). Even a Gaussian noise of root-mean-squares ± 50 m added on the bedrock topography leads to errors of a few per cent on the flow-parameter retrieval (Table 1, case 6).

Accumulation rate is established by compilations of direct measurements on stakes and by measurements of trace elements in pits or ice cores: the resulting precision is probably within 20%. An overestimation of the accumulation rate would lead to an overestimation U/H but an underestimated value of T_b , because of the importance of vertical advection. But, as shown for cases 7–8 in Table 1, the major effect is on B₀ .

At last, we analyze the role of surface topography which is very critical because it enters into the computation of the stress (Equation (2)) and of the width and direction of flowlines (Equation (8)). The error in ice topography retrieved by satellite altimetry is due to data noise, atmospheric propagation error, orbit error and slope error. Inverse techniques allow us to take into account the whole error budget as well as the a priori covariance of the signal. This leads to the unbiased computation of a precise topography and its error map. The estimated precision for Antarctica, above 1500 m, is about 1 m. The resulting precision of the surface slope on a 20 km scale is about 10⁻⁴.

The sensitivity of the ice flowlines on the topography is very large: Equation (8) shows that the channel width l(x) contributes locally and in a cumulative manner to the estimation of the parameters. An overestimation of the channel width would lead to an underestimation of the local velocity and an overestimation of the downstream velocity. We have estimated the sensitivity of the flow parameters to topography by comparing the above results with similar ones deduced from the topography of Drewry (1983). As one can see in Table 1, case 9, the dependence on the temperature (k) is strongly diminished, while n is enhanced. As we are confident in the metric precision of our topography, we conclude that stress and U/H estimations are reliable for the first time, within the model assumptions.

From Table 1 we conclude that within the model assumptions the results of the preceding section are robust. For region Β, n varies between 0.9 and 1.5, k between 0.08 and 0.1. Between the various flowlines, the scatter is comparable: n varies between 0.65 and 1.25 and k between 0.1 and 0.12, and then the activation energy Q varies between 65 and 70 kJ mol⁻¹

Conclusions

From the very precise surface topography of the Antarctic ice sheet deduced from inversion of altimeter data, we determine ice-flow direction, balance velocity and basal shear stress. We focus on regions where bottom temperature is less than −10°C, and altitude between 1500 and 3000 m. From the four selected regions we show that ice-sheet profile can be reproduced rather well by setting the Glen parameter n at 1 ± 0.25 and the activation energy around 70 ± 10 kJ mol⁻¹, for a geothermal flux of 50 mW m⁻², and when normal stress is neglected. If the temperature-dependent term was omitted, one would find a higher value of about 3 for n. A sensitivity study demonstrates that the improved precision of the surface topography is essential in this approach. When applied to the classical Drewry topography, the same calculations result in a larger value for n (n = 1.8) correlated with a lower value for activation energy.

Acknowledgements

We thank Dr D. Jenssen of the Meteorology Department, University of Melbourne, for the geophysical data sets. We are grateful to Dr P. Duval of the LGGE and Dr Ν. Young of the Australian Antarctic Division for helpful discussions and comments, and to Dr M. Monnereau for his help in the numerical resolution of the thermodynamic equation.