Flow Laws For Glacier Ice: Comparison Of Numerical Predictions And Field Measurements

Abstract

Ice flow along the 20 km long strain network up-stream of the Dye 3 bore hole in Greenland is studied in detail. By solving the force-balance equations and using selected flow laws, stresses and strain-rates are calculated throughout the section of the ice sheet. The validity of the results is evaluated by comparison with the velocity profile derived from bore-hole-tilting measurements, and with observed surface strain-rates. A number of constitutive relations are tried and most predict a velocity profile at the bore-hole site that is in good agreement with that observed, if appropriate enhancement factors are used. However, there are major discrepancies between modelled and measured surface strain-rates. Use of Nye’s generalization of Glen’s flow law, or an anisotropic constitutive relation, requires unrealistically large along-flow variations in the enhancement factor. Inclusion of normal stress effects can lead to much better agreement, but it is possible that other processes, such as dynamic recrystallization or primary creep, should be included in the constitutive relation of polar ice.

Introduction

Fundamental to understanding the dynamics of ice sheets is the constitutive relation describing the mechanical response of ice to applied stress. The most commonly used relation is Nye’s generalization of Glen’s flow law (Nye, 1953; Glen, 1958) or variations on it (e.g. Hutter, 1983). These relations are supported by laboratory and field experiments, but the experiments often invoke questionable assumptions or procedures, or they provide results that are not conclusive. It appears that the constitutive relation for polar ice is not very well established.

The ideal test for a constitutive relation would involve naturally deforming glacial ice under known stress and strain-rates. Such a test, however, has been difficult to perform because, in most cases, the deformation rate can be measured only along the surface of the deforming ice mass, or in isolated bore holes, and because stresses cannot be measured. The flow of ice leading to the Dye 3 core hole, however, offers special opportunities for testing suggested flow laws.

The data collected in the Dye 3 area in southern Greenland probably constitute the most complete set that can be used for a careful study of various constitutive relations. Of critical importance is the fact that the glacier there is frozen to the bed at a temperature of about −13°C (Gundestrup and Hansen, 1984), which is less than the temperature at which water-like phases may influence the creep of ice (Hooke, 1981). This means that true sub-solidus creep is applicable. It also means that a no-slip condition can be applied at the glacier bed. Furthermore, a core has been recovered through the full ice thickness and various physical properties such as crystal size, ice fabric, and impurity content have been measured along the length of the core (Hammer and others, 1985; Herron and others, 1985). The bore hole itself has been studied for tilting and closure rates (Gundestrup and Hansen, 1984). In addition, there are three sets of ice-sounding radar surveys (Overgaard and Gundestrup, 1985; Jezek and others, 1985a; Fisher and others, 1989), as well as measurements of surface slopes and strain-rates at the surface (Whillans and others, 1984).

Taken together, these data allow us to apply the newly refined force-budget technique (Van der Veen, 1989; Van der Veen and Whillans, 1989a) to determine how well different suggested constitutive relations apply to the flow near Dye 3. This is the subject of the present contribution.

The Constitutive Relation

When ice is subjected to a constant stress, the creep rate passes through three stages. Initially, it is high during the period of primary or transient creep. It decreases rapidly with time and may be followed by a period of secondary creep during which the strain-rate remains approximately constant. Later, when sufficient strain energy has accumulated, the ice-fabric pattern alters from a random distribution of c-axes to one in which the c-axes have developed a preferred orientation. This stage of tertiary creep is characterized by increased strain-rates.

In most applications, of the three creep stages, secondary and tertiary creep are considered to be most important for glaciers and ice sheets. Transient creep does not occur within the main body of glacier ice because this ice has been deforming within a given stress regime for a sufficiently long time to reach the stage of secondary creep. Paterson (1981) suggested that the tertiary creep follows a similar relation to secondary creep but with rate factors corresponding to much softer ice. According to this suggestion, secondary and tertiary creep are alike and, following most authors, we do not make a distinction here, except to consider the rate factor as poorly constrained and to allow for anisotropic effects.

Secondary and tertiary creep are commonly described by Nye’s generalization of Glen’s empirical representation,

, where

. is the effective strain-rate, and τ_e is the effective stress (the second invariant of the strain-rate and deviatoric-stress tensors, respectively). More specifically, a particular strain-rate component (

) is linked to the corresponding deviatoric stress component (τ′ _ij as

(1)

Discussion has mainly focussed on the values of the rate factor A and of the exponent n (e.g. Hooke, 1981; Budd and Jacka, 1989).

It is widely accepted that the rate factor A is an exponential function of temperature. For temperatures less than −10°C, this dependency follows an Arrhenius-type relation A = A ₀ exp[-Q/RT], where Q is the activation energy for creep, R is the gas constant, and T is the temperature in Kelvin. Although this relation is inappropriate at temperatures above −10°C (probably because of formation of a fluid phase on grain boundaries (Hooke, 1981)), it is commonly used for higher temperatures as well, but with a different value for the activation energy Q (Paterson, 1981). The rate factor also depends on the hydrostatic pressure, although this effect is small even for pressures that exist at the base of thick ice sheets (Rigsby, 1958).

The pre-exponential factor A_Q is not a constant but depends on crystal structure, concentrations of impurities, and perhaps other factors (Hooke, 1981; Budd and Jacka, 1989). Also, effects of crystal orientation are commonly incorporated into A ₀, for instance, by introducing an enhancement factor (Shoji and Langway, 1984; Dahl-Jensen, 1985). This point will be discussed in more detail below.

Many laboratory experiments on the deformation of ice have been carried out, and the results generally support values of n around 3 for effective stresses less than 100 kPa, and larger values for stresses in excess of 1000 kPa (Hooke, 1981). This last result is not very relevant for glaciers, where such large stresses are rare. The low-stress experiments, however, must be interpreted with caution because in many of them the ice is deformed, for a limited time, to strains of only a few per cent. In reality, ice in glaciers and ice sheets has been deformed for thousands of years, to much larger strains. It has not been practical to conduct laboratory experiments for sufficient time for the ice samples to reach strains that are relevant for glacier ice.

Field measurements have also been used. The first such study is that of Nye (1953), who measured the closure rate of tunnels in the Jungfraufirn. The tilting of a near-vertical bore hole has been used many times (e.g. Nye, 1957; Shreve and Sharp, 1970; Dahl-Jensen, 1985; Doake and Wolff, 1985; Pimienta and Duval, 1987), and measurements of the spreading rate of floating ice shelves have also been used to infer constitutive properties (e.g. Thomas, 1973; Jezek and others, 1985b). A summary of results obtained prior to 1981 can be found in Hooke (1981); a more recent overview has been given by Budd and Jacka (1989).

None of the interpretations of field data, however, leads to an unequivocal constitutive relation. A major problem is that field measurements yield only strain-rates. So, in order to complete the analysis, the stresses must be inferred indirectly, and correlated with the measured strain-rates. In most cases, the stresses are calculated from a theory that involves simplifications to make the problem solvable. These simplifications, however, may be very important to the interpretation of the field data.

As an example, consider two recent interpretations of bore-hole data (Doake and Wolff, 1985; Pimienta and Duval, 1987). In these studies, the shear stress at depth is calculated from the approximate solution derived by Nye (1957), while longitudinal stresses are either neglected or calculated in an ad hoc fashion, for example, by assuming that longitudinal stretching is constant with depth. By considering logarithmic plots of the shear strain-rate

versus shear stress τ _xz , Pimienta and Duval (1987) concluded that the tilting measurements of the Dye 3 bore hole suggest a quasi-linear constitutive relation. Doake and Wolff (1985) were more cautious and concluded from data of four bore holes that linearity might be as appropriate as the commonly used non-linear constitutive relation. There are thus major ambiguities in the determination of flow parameters from bore-hole tilt rates.

The major problem with the interpretation of bore-hole tilt data is that, in practice, the effective shear stress, τ_e, various little with depth because of the role of longitudinal deviatoric stresses. As inspection of Equation (1) shows, a range of τ_e is needed to determine precisely the value of the exponent n. However, along the Byrd Station Strain Network (West Antarctica), calculations using Equation (1) show that the square of the effective stress varies typically by a factor of only 3 from the surface to the bed (Van der Veen and Whillans, 1989b). This variation is much smaller than that of the temperature-dependence of the rate factor A (about a factor of 20), which makes it difficult to design a test using bore-hole data to discriminate clearly the value of n from uncertainties in the activation energy. This is a common problem, and similar objection against the interpretation of bore-hole tilting data has been raised by Dahl-Jensen (1989), based on her modelling results for flow along flow lines.

An innovative technique to interpret the spreading of the Ross Ice Shelf was used by Jezek and others (1985b). The principal tensile stress was estimated from the height of bottom crevasses and compared to strain-rates measured in nearby crevasse-free ice. Disappointingly, however, the results show such scatter that no particular constitutive relation can be supported.

The experiments using laboratory and field measurements have, at best, determined values of the flow parameters A and n, if Equation (1) is adopted. None of the experiments has actually confirmed the validity of the general form of the relation. In fact, it can be argued that many laboratory experiments have disproved the applicability of Glen’s law. Even in 1958, Glen recognized that the few experimental results then available (most notably those from Steinemann) disagreed with the simple flow law commonly used in glaciology. Therefore, he proposed a more general relation, namely (Glen, 1958)

(2)

where IΔ₂ and IΔ₃ are the second and third invariants of the deviatoric stress tensor, and δ _ij is the Kronecker delta. Because this expression contains two unknown stress functions (F ₁ and F ₂), uniaxial stress experiments (as most laboratory experiments are) are insufficient to determine both functions. As discussed by Morland (1979), combined shear and compression experiments are required to resolve the form of the two unknown stress functions. Unfortunately, this point has been taken up by very few experimentalists (cf. Baker, 1987). Rather, the usual procedure is to neglect the second term on the right-hand side, and also to ignore the dependency of the first term on the third invariant of the stress deviator.

There is a further fundamental question concerning the validity of applying Equation (1) to natural ice as found in glaciers. The fabric pattern changes gradually from a random distribution of the c-axes near the surface to a pattern with c-axes clustered around a - few directions at greater depths (Kamb, 1959; Rigsby, 1960; Gow and Williamson, 1976; Herron and others, 1985). This change in fabric is caused by rotation of basal planes associated with deformation, growth of certain ice grains, and recrystallization (Matsuda and Wakahama, 1978; Hooke and Hudleston, 1980; Azuma and Higashi, 1985; Alley, 1988). As a result, the mechanical properties of polar ice can be expected to be a function of its deformation history (Alley, 1988).

Laboratory tests indicate that the rate of deformation for a given shear stress depends strongly on the orientation of the c-axes. Experiments carried out on the Wisconsin ice retrieved from Dye 3 show that shear strain-rates, τ _xz , in this ice, with a pronounced clustering of c-axes along the vertical z-axis, are an order of magnitude larger than in randomly oriented laboratory ice at the same stress level (Shoji and Langway, 1984, 1985).

The usual procedure to incorporate these effects of anisotropy into the constitutive relation is by introducing an enhancement factor. To this end, the rate factor A is written as the product of a temperature-dependent part A(T), as discussed above, and the depth-dependent enhancement factor E(z): A = A(T)E(z) (Dahl-Jensen, 1985). This is an expedient technique, but accounting for anisotropic effects through adjustment of the rate factor is not really satisfactory.

Enhancement factors do not account for anisotropy. In ice with c-axes aligned with the vertical axis, it is reasonable to expect that the ice will be stiffer to longitudinal stresses, and softer with respect to shear stresses, than isotropic ice. Such a distinction, however, cannot be described by Equation (1) with a single depth-varying rate factor. Of course, this criticism can be countered by introducing different rate factors for different strain-rate components.

A theoretically more sound constitutive relation that includes anisotropy has been formulated by Johnson (1977), and was, in modified form, applied to glacier ice by Lliboutry and Andermann (1982) and Pimienta and others (1987). With some minor modifications, these works are followed here, in an effort to determine a more satisfactory constitutive relation for polar ice.

Other forms of anisotropic constitutive relations for polar ice have been proposed (Lile, 1978; Baker, 1981, 1982; Jacka and Budd, 1989). Most of these were derived from laboratory experiments and are not yet fully supported by a mathematical theory. Rather, we start out with Johnson’s most general form of the anisotropic constitutive relation and study its effect on glacier flow.

A full two-dimensional force-balance study is conducted. Stresses and strain-rates along a 20 km stretch of flow line up-stream of the Dye 3 bore hole (Greenland) are calculated results, and are expected to accord with data from along the surface and with bore-hole data.

Force Budget: Method Of Solution

The ice flow along a small grid up-stream of the Dye 3 bore hole is considered. The methods described in Van der Veen (1989) and Van der Veen and Whillans (1989a) are used to calculate all relevant stresses and velocities at depth. In short, calculations start at either the bed or top surface and proceed through the glacier in small increments until the other surface is reached.

For the first calculations, an upward scheme is used (Van der Veen, 1989). An initial estimate for the basal drag, and prescribed basal velocities (zero in this case) are used to calculate vertical shearing at the base, from which velocities at some depth above the bed are calculated. These velocities are used to compute the forces on vertical faces of a column extending from the bed to that depth. From force balance, the shear stress at this depth can then be calculated. This vertical shearing is used to compute the velocities at the next higher layer. Thus, by progressively going upward, velocities and stresses are calculated throughout a section of the glacier. When the surface is reached, the net stress at the surface may not be zero, as it physically should be. If it is non-zero, an improved value for basal drag is computed, and the upward calculation is repeated. This iteration cycle is carried out until the calculated surface stress is near zero (less than, say, 1 kPa).

In later calculations in which the properties of the basal soft-ice layer are varied along flow, the force-budget technique is reversed and calculations proceed from the top down (Van der Veen and Whillans, 1989a). This method is more direct and faster but does not objectively predict zero basal velocity, as must happen in Nature. Rather, the necessary enhancement factor for deep ice is determined in order to satisfy zero basal velocity. Neither solution scheme involves arbitrary assumptions other than that of plane flow. This assumption is necessary because the transverse dimension of the strain grid at Dye 3 (4 km with three measurement sites) is too small to warrant full three-dimensional calculations. It also seems appropriate because comparison of surface velocities along the three columns of stakes in the strain grid shows that the flow is almost parallel (Whillans and others, 1984), so the plane-flow assumption is probably not unrealistic.

In the first, upward, calculations, two important results are the basal drag and the surface velocity along the length of the strain net. The velocities can be compared with those measured (Whillans and others, 1984), as a test of the results of the calculations.

Data

The Dye 3 grid extends 20 km along the mean direction of ice flow, and consists of three parallel rows of poles, about 2 km (one ice thickness) apart (Whillans and others, 1984). The bore hole is located at the down-flow end of the grid (Fig. 1).

Surface elevations were measured by barometric altimetry, referred to the two satellite-tracking stations (Whillans and others, 1984). Elevations are believed to be least accurate at the up-stream end of the grid, where the error may be as large as ±2 m. Measured surface elevations were slightly smoothed using a Gaussian smoothing scheme to ensure that there are no short-scale variations that could adversely affect the solution. The smoothed elevations relative to the surface elevation at Dye 3 are shown in Figure 2.

Fig. 1. Map of Greenland showing the location of Dye 3 strain net, and detailed map of the Dye 3 area (after Dahl-Jensen (1985), based on data from Overgaard and Gundestrup (1985)). The three straight lines represent the three grid lines and the dashed line the estimated flow line leading to the Dye 3 bore hole (cf. Whillans and others, 1987). Cross-lines indicate stations along the central line of the strain grid used in this study. Open circles denote satellite-tracking stations (Whillans and others, 1984).

Fig. 2. Measured surface elevations relative to the elevation at the bore hole (from Whillans and others, 1984) and basal topography (from Jezek and others, 1985a; Overgaard and Gundestrup, 1985) along the center line of the strain grid.

Basal elevations, also shown in Figure 2, were letermined by radar sounding (Jezek and others, 1985a; Dvergaard and Gundestrup, 1985). Although the accuracy of;ach radar survey is about ±15 m, the relative difference between bed elevations from the two data sets is as large as 100 m at a few places. These discrepancies have no;ignificant effect on the results presented here.

Because the rate factor in the constitutive relation is :emperature-dependent, the vertical temperature profile is ilso needed. For this, the profile measured in the bore hole Gundestrup and Hansen, 1984) is used (Fig. 3). This profile is applied to the entire grid. Considering the short length of the flow-line segment (20 km), and the fact that the results of the force-budget calculations are rather insensitive to the value of the temperature at depth (Van der Veen and Whillans, 1989b), this approximation is believed to have only minor effects on the results discussed below.

Fig. 3. Temperature measured in the Dye 3 bore hole (from Gundestrup and Hansen, 1984) and enhancement factors at depth calculated by Dahl-Jensen (1985) for the Dye 3 bore hole. The vertical scale gives the dimensionless depth, s = (h - z)/H, in which h represents surface elevation, H ice thickness, and z elevation.

A small problem arises because the detailed strain grid leads up to, but does not enclose, the bore-hole site, so accurate surface elevations are not available beyond the bore hole. However, for reliable force-budget calculations at the bore-hole site, the numerical grid should extend beyond the bore hole. The more widely spaced surface elevations given in Whillans and others (1987) do enclose the bore hole and were used to extrapolate the surface slope. Basal elevations were also extrapolated but these have little effect on the outcome of the force-budget calculations. In the following, only results that pertain to the actual strain grid are shown.

This extrapolation of the surface could lead to slightly inaccurate force-budget calculations at the edges of the model domain, so that the velocity-depth profile computed for the bore-hole site is not entirely correct. However, the main problems associated with matching model results with observations are not associated with the observed bore-hole tilting, but rather with the surface variations in velocity.

Depth Variation In Ice Properties

Analyses of ice samples from the bore hole have shown that the physical and chemical properties vary considerably with depth. On the basis of these properties, the entire column can be divided into three layers with marked differences (Hammer and others, 1985; Herron and others, 1985).

The bottom 25 m (1.2% of the ice thickness) consists of silty ice with a high dust concentration and small crystal diameters. On top of this basal layer lies the 230 m (11.3%) thick layer of Wisconsin ice. This ice features high dust concentrations, small crystal diameters, and a strong single-maximum vertical orientation of the c-axes. In the zone of Holocene ice (the upper 1786 m, or the top 87.5%), the concentration of dust is an order of magnitude smaller than in the underlying ice, and the ice fabric changes gradually to a random orientation of the c-axes at shallow depths.

Mechanical tests indicate that the Wisconsin ice is much softer than the Holocene ice (Shoji and Langway, 1984, 1985, 1988). This means that, for a given shear stress, the rate of deformation is larger in the Wisconsin ice than in the Holocene ice. To account for this effect, it is customary to include an enhancement factor as discussed above (Dahl-Jensen, 1985; Dahl-Jensen and Gundestrup, 1987).

As argued above, accounting for anisotropy effects through an enhancement factor is, from a theoretical point of view, not very attractive. Therefore, the constitutive relation formulated by Johnson (1977) is also applied to the Dye 3 strain grid. In Appendix A, a derivation of this constitutive relation is given. It is the simplest generalization from Equation (1) for steady creep in polycrystalline ice to describe creep in ice with a pronounced preferred orientation of c-axes. Two material parameters (α and β) account for the constitutive response to, in this case, horizontal normal stresses. A third parameter (γ) describes the response to shear perpendicular to the c-axes. Values of these parameters have not been determined for polar ice, but γ may be taken as analogous to the enhancement factor of Dahl-Jensen (1985).

The anisotropic constitutive relation derived in Appendix A is only applicable to glacier ice in which the c-axes are aligned with the vertical z-direction. A more general theory should also include an equation describing the change with time of the direction of anisotropy. However, for glacier ice, such a theory has not been developed. As argued in Appendix A, a reasonable first step is to investigate theories for anisotropy that can explain the observed pattern of flow for the current situation in which the preferred c-axis orientation is along the vertical.

Modelling The Ice Flow

The simplest model of a glacier is the laminar or lamellar model (e.g. Paterson, 1981, p. 85–89) in which a balance between driving stress and basal drag is assumed. Also, in this model, vertical variation in horizontal velocity depends on the shear stress, τ _xz , only (with x horizontal in the direction of flow, and z vertically upwards). Because the driving stress varies linearly from zero at the surface to the maximum value at the bed, the velocity can be expected to be large where surface slopes are large, and vice versa. As shown in the upper panel of Figure 4, the driving stress shows very large short-scale variations along the strain grid. Consequently, velocities calculated from the laminar-flow theory would display enormous, and unrealistic, horizontal variations.

Fig. 4. Driving stress for the full thickness and calculated basal drag along the Dye 3 flow line, using Glen’s flow law (n = 3, A₀ = 4.0 × 10⁵ kPa⁻³ year⁻¹, Q = 60 kJ mol⁻¹) including the enhancement factors shown in Figure 3, and using the anisotropy law of Johnson (1977).

A more realistic flow model is one that includes the effects of gradients in longitudinal stresses on the shear stress at depth. Part of the variation in driving stress is balanced by horizontal stress gradients, and the result is that basal drag varies much less than does the driving stress. As a result, surface velocities vary more slowly than according to the laminar model. This is a major improvement.

This model does not take into the account the markedly different constitutive properties of the softer basal Wisconsin ice. This can be treated through the use of an enhancement factor for deep ice by allowing for crystal anisotropy by varying properties along flow, or by including normal stress effects.

Enhancement-factor model

This model is very similar to that used by Dahl-Jensen (1985). The force-budget technique (Van der Veen, 1989), together with the common flow law, but including an enhancement factor, is applied to calculate stresses and velocities at depth. The enhancement factor shown in Figure 3, scaled to the ice thickness, is applied along the entire flow line. The rate factor is assumed to be temperature-dependent according to the Arrhenius relation (Q = 60 kJ mol⁻¹, n = 3, and A ₀ = 4.0 × 10⁵ kPa⁻³ year⁻¹), with the temperature profile shown in Figure 3. Results of the calculations are presented in Figures 4, 5, and 6 (curves labelled “Enh”). As expected, basal drag varies less than the driving stress (Fig. 4), because of the role played by horizontal stress gradients.

Fig. 5. Calculated velocity profile at the bore-hole site. The curve labelled “Enh” corresponds to the calculation based on Glen’s law including enhancement factors, and that labelled “An” to the calculation based on Johnson’s anisotropic law. Velocities derived from bore-hole-tilting measurements (from Gundestrup and Hansen, 1984) are indicated by crosses. Velocities have been non-dimensionalized as fu(s) [u(s) - u(0)]/u(0).

Comparison of the calculated velocity profile at the down-stream end of the grid with that derived from bore-hole tilting measurements (Gundestrup and Hansen, 1984) shows that agreement is satisfactory. This is to be expected because the enhancement factor used in this calculation is derived from a similar calculation by Dahl-Jensen (1985), who matched the calculated profile with that observed.

Although this model of full force budget with enhancement factor is a major improvement over the laminar-flow theory, it is not sufficient. As shown in Figure 6, the calculated variation in surface velocity is much larger than observed. Therefore, a further level of sophistication is required. Dahl-Jensen encountered the same problem with her calculations, but did not emphasize the disagreement and its consequences, as we do here.

Fig. 6. Surface stretching and surface velocities along the Dye 3 flow line. The heavy curves are the results of force-balance calculations, while the thin curves are data from Whillans and others (1984).

Anisotropy

Accounting for the softer Wisconsin ice through an enhancement factor is insufficient to bring calculations in line with Nature. A higher level of sophistication is to use a mathematically sound constitutive relation that describes the anisotropic properties of, especially, the deeper Wisconsin ice. Such a constitutive relation has been derived by Johnson (1977), and that work is followed here with some minor modifications (cf. Appendix A).

The curves labelled “An” in Figures 4, 5, and 6 show the results of a force-budget calculation based on this anisotropic constitutive relation. Although an enhancement factor is not used explicitly, it is effectively replaced by the parameter γ that describes the stiffness or softness to shear of one layer over another. It has the same pattern of depth variation as shown in Figure 3 for the enhancement factor used in the previous model. The other two parameters appearing in Johnson’s law (α and β) are taken to be equal to their isotropic value, 1.

Inspection of the calculated velocity profile at the bore-hole site (Fig. 5) shows that the horizontal velocity has a small minimum at about 80% depth (just above the Wisconsin-Holocene ice transition). This is very likely an artefact of the model caused by the rather sharp transition boundary, and no great significance should be attached to it. The match in this figure could be greatly improved by using a more complex vertical variation of the anisotropy factor, γ. Model calculations using different depth variations of γ indicate such a change would have little effect on the calculated patterns of basal drag and surface stretching in Figures 4 and 6. Therefore, no attempt was made to improve the agreement between the calculated and measured velocity profile at the bore-hole site.

Many experiments were conducted with different degrees of anisotropy, but the results presented in Figures 4 and 6 are about the best match with surface observations that was obtained (using γ = 1.2 × E at depths greater than 85% of the local ice thickness, and γ = 2 (the isotropic value) at other depths; both α and β are kept fixed at their isotropic value, 1). Comparison of the curves labelled “An” with those labelled “Enh” in these figures shows that the anisotropic model yields virtually the same results as the enhancement model. The role of longitudinal stress gradients in the balance of forces is still too small, and calculated variations in the surface strain-rate pattern are too large compared to those measured.

Accounting for the soft Wisconsin ice through an anisotropic-based constitutive relation thus fails to model the ice sheet correctly. To obtain more satisfactory results, either a more complicated constitutive relation must be chosen, or the properties of the soft-ice layer must be allowed to vary along the flow line. In the experiments discussed so far, the enhancement factor, or the anisotropy factor, is assumed to be laterally uniform. That is, the same factors were, scaled to the ice thickness, applied to the entire segment of the flow line. However, the Wisconsin ice may vary in the direction of flow, and this possibility leads to the next model.

Along-flow variations

By horizontally varying the thickness of the Wisconsinan layer, the enhancement-factor, or the anisotropy to shear, agreement between calculated and measured surface velocities and strain-rates can be improved.

There is, however, a practical difficulty in using the “upward calculation” scheme for such calculations. The introduction of horizontal gradients in ice properties affects horizontal gradients in stress (and hence basal drag) in a complex way, so that it becomes difficult to anticipate how the enhancement factor, or the thickness of the soft-ice layer must vary. The downward calculation scheme described in Van der Veen and Whillans (1989a) does not pose these problems, and is therefore used instead.

The basic principle of the downward scheme is the same as that of the upward scheme, except that calculations start at the surface instead of at the base. At the surface, the measured velocities are prescribed and used to calculate shear strain-rates at the surface. These are used to compute velocities at some depth layer below the surface and, together with the force-balance equation for that level, this allows shear strain-rates to be calculated. Calculations proceed gradually downward until the top of the Wisconsin ice layer is reached. In this layer, the horizontal ice velocity varies in a prescribed fashion (linearly) from the calculated value at the top of the soft-ice layer to zero at the bed. This assumption immediately yields all strain-rates as well as the vertical component of velocity. The enhancement factors within the Wisconsin ice required to satisfy force balance can then be calculated. Alternatively, the depth to the Holocene-Wisconsin transition could be varied for constant enhancement factors, but in the numerical model this is not very convenient.

The downward calculation scheme was tested extensively, primarily using data from the Byrd Station Strain Network (Van der Veen and Whillans, 1989b). One important issue, not discussed in that paper, concerns the numerical stability of the scheme. Tests described in Appendix C show that no numerically induced variations or instabilities develop during the downward calculation. Thus, the calculated pattern of basal drag is real.

Figure 7 shows the results of a calculation based on Nye’s generalization of Glen’s flow law, with the depth of the Wisconsinan layer kept fixed at 15% of the ice thickness along the entire flow line. Because the more slowly varying measured longitudinal strain-rate gradients are used, the calculated variations in basal drag are much larger than in Figure 4.

Fig. 7. Results of a downward force-balance calculation using enhancement factors, with the depth of the Wisconsin ice layer fixed at 15% of the ice thickness along the entire flow line, and prescribed zero basal velocities. The upper panel shows calculated basal drag, the middle panel enhancement factors required to satisfy force balance, and the lower panel the basal elevation with the linear trend subtracted (n = 3, A₀ = 4.0 × 10⁵ kPa⁻³ year⁻¹, Q = kJ mol⁻¹).

The calculated enhancement factor shows variation over two orders of magnitude, along the flow line. Comparison of the middle curve in Figure 7 with the curve of the driving stress (Fig. 4) shows that the enhancement factors are large where the driving stress is small, and vice versa, as expected in order to obtain the observed slowly varying surface velocities.

The lower curve in Figure 7 shows the basal elevation with the linear trend removed. It shows a strong anti-correlation with the curve above. Where the bed is deep, enhancement factors are large whereas over basal highs calculated enhancement factors are small. If, instead of calculating horizontally varying enhancement factors, the thickness of the soft layer had been varied, the corresponding result would be thick, soft ice in the basal lows. Thus, basal lows appear to contain especially soft ice, either because the thickness of the soft-ice layer is larger or because the basal ice is softer, or a combination of both.

A disturbing feature in the results shown in Figure 7 is that the variation in the enhancement factor is very large. In fact, when interpreted in terms of thickness of the soft-ice layer, the two peaks at 11 and 18 km imply a soft-ice layer extending almost to the surface (using the value of the enhancement factor obtained at the bore hole; Fig. 3). This is not realistic and there appears to be a major problem with the assumptions on which this calculation scheme is based.

Using the calculation scheme based on the anisotropic constitutive relation in a similar fashion to calculate the parameter γ, describing enhanced deformation under shear, gives virtually the same results. As shown in Figure 8, large variations in γ are required to satisfy force balance under the constraints imposed at the upper and lower boundaries.

In order to obtain more realistic results (that is, more slowly varying enhancement factor or anisotropy factor), the role of the longitudinal gradient terms must be enhanced. In the “enhancement” model, this can be achieved by decreasing the value of the rate factor A ₀ (stiffer ice). However, experiments conducted show that this has only a minor effect. In the “anisotropy” model, longitudinal gradient terms become more important when α and β are adjusted to make the ice stiffer (smaller values of α and β) than in the isotropic case.

Results of two calculations in which α increases linearly from a minimum at the surface to the isotropic value at the top of the Wisconsin layer, and remains constant below (β = α), are also shown in Figure 8. The effect is considerable and variations required in the shear parameter γ are much less than before. In fact, γ varies only by about a factor of 4, which seems to be entirely within reason.

Although the results depicted in Figure 8 suggest that the “anisotropy” model, including along-flow variations of the Wisconsin layer, is a realistic model for the ice flow in the Dye 3 area, a worrying issue is that the value of α in the upper ice layers is extremely small. To obtain realistic results, the upper ice must be more than 500 times stiffer than commonly assumed. Measured differences in softness between the Holocene and the Wisconsin ice are typically only about a factor of 5–10 (Shoji and Langway, 1988). Furthermore, optical studies of c-axis orientation show that the upper layers are random and therefore are isotropic (Herron and others, 1985). The results shown in Figure 8 do not have a strong physical basis, and appear to be fortuitous at best.

Previous authors also encountered problems in modelling the constitutive properties of the upper ice layers (Shreve and Sharp, 1970; Raymond, 1973; Hooke and Hanson, 1986). Raymond (1973) concluded that the upper one- to two-thirds of Athabasca Glacier constitutes an anomalous zone in which there is either a strong effect from a complex distribution of stress arising from the longitudinal stress gradients, or more complicated rheological behavior. Similarly, Hooke and Hanson (1986) indicated that Equation (1) may not be applicable near the glacier surface where a low-stress situation prevails in which no single stress component is clearly dominant.

It seems that an additional process that is important on a short scale (less than about two ice thicknesses) needs to be included in the model for ice flow. It is not clear what this mechanism may be, but one possibility is that of normal stress effects. These effects are known to be important in other disciplines that consider non-Newtonian fluids (e.g. Schowalter, 1978).

Normal stress effects in polar ice

In 1945, Reiner predicted on purely mathematical grounds that normal stress effects (which he referred to as dilatancy, a term coined by Reynolds in 1885) should be present in any viscous material whose rheological properties are described by a non-linear functional relationship between deviatoric stress and strain-rate. Two years later, Weissenberg (1947) reported observations of such effects in a variety of materials (saponified oils, solutions of rubber, starch, cellulose acetate, etc.) that were subjected to simple shear. In glaciology, this phenomenon has received very little attention, although it was mentioned by Glen (1958).

As noted earlier, the most general form of the constitutive relation is (Reiner, 1945; Glen, 1958)

(2)

where F ₁ and F ₂ are unknown functions of the second (I′₂) and third (I′₃) invariants of the deviatoric stress tensor. The second term on the right-hand side (commonly neglected in glaciology) is responsible for the manifestation of normal stress effects. Consider, for example, the situation in which a simple shear stress τ _xz is applied. From Equation (2), the strain-rate tensor is found to be

(3)

The second term indicates a swelling or contraction perpendicular to the plane of shear, given by

(4)

and corresponding contraction or swelling in the plane of shear, so that the ice volume does not change. This expression suggests that the ice will attempt to expand or contract in the presence of shear, but appropriate deviatoric normal stresses may develop to oppose this tendency.

Fig. 8. Basal drag and shear parameter γ required to satisfy force balance, when zero velocities are prescribed, and the depth of the Wisconsin layer is fixed at 15% of the ice thickness, along the entire flow line. All three curves are the result of a downward calculation based on Johnson’s anisotropy law. The longitudinal stretching parameter a was assigned the isotropic value 1 (curves labelled 1), or increased linearly from a minimum value at the surface (0.05 for curves 2, and 0.005 for curves 3), to the isotropic value at the top of the Wisconsin layer (and constant below). The value of the rate factor used in all three calculations is A₀ = 4 × 10⁵ kPa⁻³ year⁻¹.

Two papers that discuss models incorporating normal stress effects into glaciology are those of McTigue and others (1985) and Man and Sun (1987). The former study adopted as a tentative model for polar ice the second-order fluid. The proposed constitutive relation, as noted by the authors, cannot be simplified to Nye’s generalization of Glen’s flow law (Equation (1)). Therefore, the two modifications to the usual flow law (the modified second-order fluid, and the power-law fluid of grade 2) introduced by Man and Sun (1987) are used here. Appendix B presents these two constitutive relations, and their application to plane flow.

The anisotropic law developed in Appendix A also includes normal stress effects, for the case in which α is not equal to β (see Equation (A21) in Appendix A). However, to keep the different physical processes clearly distinct, it seems preferable to use the relations given by Man and Sun (1987), rather than to include normal stress effects resulting from anisotropic behavior.

McTigue and others (1985) asserted normal stress effects may affect the shape of the free surface of a glacier in an open channel, and may account for a depression of more than 1 m in the surface. Man and Sun (1987) found that the free-surface depression (or heave) is even smaller than this. Despite the apparent smallness of these effects, it is worthwhile investigating the applicability of Man and Sun’s models to the flow of polar ice.

Fig. 9. Enhancement factors required to satisfy force balance, when zero basal velocities are prescribed, and the depth of the Wisconsin layer is fixed at 15% of the ice thickness, along the entire flow line. The upper curve applies to a downward calculation based on Nye’s generalization of Glen’s flow law (n = 3, A₀ = 4.0 × 10⁵ kPa⁻³ year⁻¹, Q = 60 kJ mol⁻¹), the middle curve to that using flow law I of Man and Sun (1987) (A₁ = 100), and the lower curve pertains to the calculation based on flow law II of Man and Sun (1987) (A₁ = 1.15).

Normal stress effects alone cannot account for the velocity pattern measured along the grid, and so the enhancement factor is also invoked. As in the preceding model, the depth of the Wisconsin-Holocene transition is kept fixed at 85% of the ice thickness. Zero basal velocities are prescribed and the enhancement factor required to satisfy force balance is calculated.

Figure 9 shows the calculated enhancement factor (average value for the Wisconsinan layer) along the Dye 3 flow line. The upper curve is the same as in Figure 7. Using flow law I of Man and Sun (1987) (modified second-order fluid), the curve becomes much smoother, but still with one very large peak, about 9 km up-stream of the bore hole. This peak reduces significantly when the second flow law (power-law fluid of grade 2) is used. Instead of extreme values differing by a factor of 150 (upper curve), the range has been narrowed to a factor of 20.

Further experiments are possible. For example, a different depth variation of the second rate factor (α₁) in Man and Sun’s constitutive relations could lead to less longitudinal variation in enhancement factor. However, there are no good, independent estimates of the magnitude of normal stress effects on glaciers, or of the spatial variation in thickness or enhancement factor of the Wisconsin ice layer. The range in calculated enhancement factors could be narrowed further through selection of appropriate normal stress-effect parameters but, without data against which to test the relations, there seems little point in trying to obtain “better” results.

Other possible models

There are other processes involved in creep that are not tested here. Among these are memory effects associated with the development of ice fabric and texture, processes associated with the third invariant of the stress or strain-rate tensor, and dynamic recrystallization.

The development of fabrics in ice sheets depends on the deformational history. This means that local stresses are related not only to local strain-rates, but to previous strain-rates as well. So, it may be argued that the constitutive relation should include a memory function. This would entail combining a model like that of Alley (1988) to calculate the development of fabric, together with the model used here, to solve the force-balance equations.

In order to achieve agreement between calculated and measured velocities, the memory function would have to make the ice alternately stiff and soft with distance along the flow line. The pattern should be similar to that of the calculated enhancement factor shown in Figure 7. However, all deep cores that have been studied so far indicate rather similar depth variation of fabric, and it seems unlikely that this is fortuitous. Thus, the available evidence does not favor the large longitudinal fabric variations necessary to explain the measured surface strain-rates.

Another effect could be the role played by the third invariant of the stress or strain-rate tensor. The experiments of Baker (1987) suggest that strain-rates may be dependent on the third invariant of the deviatoric stress. However, it is unclear how this effect should be implemented in the constitutive relation.

A further process is dynamic recrystallization. Jonas and Müller (1969) found this to have a large effect. When ice undergoes dynamic recrystallization, the grains or crystals that are being deformed can be continuously replaced by new and undeformed crystals. This means that the constitutive relation should take into account the softening of the migrating grain boundaries and the replacement of the deformed grains by strain-free ones, and of the rapid reduction in dislocation density. Jonas and Müller (1969) concluded that dynamic recrystallization can occur after a critical strain (depending on the stress and temperature) has been reached. At intermediate stresses, the recrystallization is periodic, leading to periodic increases in strain-rates while, under higher stresses, recrystallization can become continuous, leading to a single increase in strain-rate. If this is an important process in glacier flow, it may very well be that this, together with primary creep of newly formed crystals, should be included in the constitutive relation.

Another concern is that of three-dimensional effects at depth. The three ranks of stations constituting the strain grid are insufficient to conduct a full three-dimensional force-budget calculation of stresses and velocities at depth. It has therefore been necessary to make the assumption that the direction of flow does not change with depth. This may not be realistic and there could be important sideways deviations of the flow near the glacier bed.

Complex flow at depth would be associated with an altered stress distribution, and with larger strain-rates. These would lead to a larger effective strain-rate, and hence to softer ice. The effect is apt to be minimal, however, because velocities must approach zero at the bed so that there is an upper limit to the magnitude of strain-rates. Also, these deep strain-rates have little effect on the force-budget calculations because of the relatively high temperature at depth, which makes the basal ice soft. Thus, flow deviations can be caused by very small values of deviatoric stress. Such small stresses acting over the limited thickness of warm ice have little effect on the full-column force budget, and thus do not impact the calculated velocity variations in an important way.

Conclusions

Force-budget calculations of the kind presented here for the Dye 3 strain net are an excellent means for testing constitutive relations. To our knowledge, the method makes no important assumptions, except, of course, about the form of the constitutive relation being tested. These tests are the most discriminatory tests that we know of for the constitutive relation of naturally flowing glacier ice.

Not all suggested constitutive relations have been tested here. That would be a major undertaking. Rather, the most commonly used relation is tested, together with some of the adjustments suggested to allow for softer Wisconsin ice, for ice-crystal-orientation fabric, normal stress effects, and possible longitudinal variations in the thickness of the soft Wisconsin-age ice. Most tests failed to provide satisfactory agreement between model calculations and observations, and it must be concluded that the standard constitutive relation, with adjustments for enhancement effects or crystal anisotropy, is deficient.

Inclusion of normal stress effects can account for much of the observed behavior. Using flow law II of Man and Sun (1987), describing the constitutive response of a power-law fluid of grade 2, the ice becomes sufficiently non-lamellar to make longitudinal stress gradients more important. This means that, in order to explain the stretching rates measured at the surface, only moderate along-flow variations in enhancement factor, or thickness of the soft Wisconsin ice layer, are required. In all other tests described here, either the calculated surface strain-rates are much larger than those observed, or unrealistically large variations in enhancement factors are needed.

The tests conducted here are for flow at the 2 km scale, and calculated results are sought that match both vertical and horizontal variations in velocity. For certain applications, such tests may be overly demanding. For example, the large-scale form of ice sheets can be reasonably well reproduced using the conventional flow law (cf. Paterson, 1981). Also, it should be emphasized that qualitative agreement is obtained here, and that the predicted strain-rate variations are only 2–3 times larger than observed. Nevertheless, the apparent failure in explaining the flow along the Dye 3 grid is important.

Recommendations can be made for future such studies. The Dye 3 strain grid is barely sufficient, and it would be better if it extended well beyond the bore-hole site, and if it contained more stations in the cross-flow direction. More importantly, it would be very helpful to have better data on internal radio-reflecting layers. The layers traced along the strain grid extend to only half of the depth (Whillans and Jezek, 1987). This depth limit is believed to be equipment-controlled because another radar obtained deeper layers at other similar sites (Gudmandsen, 1975). If deeper layers were available, they could be used as an additional test of the model, and perhaps as an independent means of testing for complex three-dimensional flow effects at depth.

Acknowledgements

Valuable criticism from J. Bolzan, R. Hooke, C. Doake, and L. Morland helped improve this manuscript. We wish to thank D. Dahl-Jensen for the hospitality shown to C.J. van der Veen while in Copenhagen, and for providing us with the values of enhancement factors. The U.S. National Science Foundation, Division of Polar Programs, is thanked for their support of this study through grant DPP-8716447 to I.M. Whillans.

This is contribution No. 721 of the Byrd Polar Research Center, The Ohio State University.

Appendix A

Tertiary Creep In Axially Symmetric Ice

Mathematical structure of Nye’s generalization of Glen’s law

To arrive at a constitutive relation describing creep in axially symmetric ice, the approach of Johnson (1977) is followed. Before including anisotropic effects, however, it is instructive to take a closer look at the mathematical structure of the conventional constitutive relation.

For isotropic ice, the constitutive relation, linking the strain-rate tensor D to the deviatoric-stress tensor T′ reads (Paterson, 1981, p. 84–85):

(A1)

with the effective stress defined as

(A2)

To cast Equation (A1) in a different form, we need some tensor relations.

Let ƒ be an arbitrary scalar function of the full-stress tensor T. Defining the tensor ∂ƒ/∂T, with components ∂ƒ/∂T _ij , it can readily be shown that the following equality holds (I being the unity tensor):

(A3)

A second useful relation is

(A4)

Taking ƒ(T) = τ_e in Equation (A3) and using Equation (A4), it follows that

(A5)

Note that, from the definition of deviatoric stresses, tr(T′) = 0, so that the second term on the right-hand side of Equation (A3) is zero.

Using relation (A5), the isotropic constitutive relation (A1) can also be written as

(A6)

with

(A7)

The scalar function Φ is usually called the dissipation potential, and it is proportional to the power dissipated per unit volume. It is important to note that this physical meaning of Φ implies that it must always be positive or zero.

Noting that the second invariant of the stress-deviator tensor is

the dissipation potential Φ can also be written as

(A8)

In other words, conventional wisdom has it that Φ is ony a function of the second invariant of the stress-deviator tensor (note that the first invariant, I ₁(T′) = tr(T′), is zero, by definition). Baker (1987) concluded from the results of his multiple-stress experiments that the third invariant I ₁(T′) = det(T′)) might also be important. However, his experimental results did not demonstrate the proportionality between the deviatoric-stress components and the corresponding strain-rates. Keeping in mind that the analysis presented below allows inclusion of this effect rather simply, it is assumed here that the third invariant of the stress-deviator tensor enters into neither the isotropic nor the anisotropic constitutive relation.

General form of the anisotropic constitutive relation

In this study, attention is restricted to transversely isotropic ice, or equivalently, ice with cylindrical symmetry. This means that with the axis of symmetry in the z-direction, say, the ice has a different stiffness with respect to shear stresses (τ _xz , and τ _yz ) than with respect to stresses in the other directions (τ _xx , τ _yy , or τ _xy ), but the stiffness is the same for both stresses τ _xx and τ _yy . Ice possessing such properties of symmetry is commonly found in polar ice sheets.

The orientation of the axis of symmetry is denoted by the unit vector n. Define the tensor

(A9)

with components

(A10)

Because, for polar ice the axis of symmetry is usually the vertical z-direction, we choose n = (0,0,1). Then the only non-zero component of M is M_zz = 1.

We assume that the anisotropic constitutive relation can be written in the general form of Equation (A6). If Q (with transpose Q ^T) is an orthogonal matrix representing either a rotation around the n-direction or a reflection in a plane containing this direction, then, for transversely isotropic materials, the dissipation potential Φ must meet the requirement

Ericksen and Rivlin (1954) have shown that this is only the case when Φ is a function of the five scalars

As noted above, the effect of the third invariant of the deviatoric-stress tensor on the constitutive relation will not be included here. Furthermore, the first invariant is zero by definition, thus Φ is a function of the following three scalars only:

(A11)

It is now convenient to introduce the following three invariants:

(A12)

These invariants are similar to those used by Lliboutry and Andermann (1982). Instead of T ₂ ² as defined above, however, they used

Contrary to what Lliboutry and Andermann (1982) claimed, their three invariants do not add up to the square of the effective stress

because of the cross-product τ′ _xx τ′ _yy in their

Of course, any combination of the three invariants can be used to construct the dissipation potential Φ. Here, we adopt the heuristic approach of choosing the simplest possible form that is reducible to expression (A7) for the isotropic case. Hence, we use

(A13)

(A14)

or, equivalently,

Again, our approach differs from that of Lliboutry and Andermann (1982), who adopted an expression that also includes cross-product terms, such as

. They argued that their expression agrees better with the few available measurements than does Equation (A14). The disadvantage of using the expression given by Lliboutry and Andermann (1982) is that it has seven unknown parameters of which only three have been measured directly while the remaining four are determined from theoretical considerations that are not clearly explained (Pimienta and others, 1987). Given the lack of observations, it seems more appropriate to use the simpler relation (A 14).

To arrive at workable expressions for the strain-rates, in addition to Equation (A3) and (A4), two more tensor relations are needed. Namely,

(A15)

and

(A16)

Using these expressions, we find

(A17)

With the axis of rotation in the z-direction, the strain-rate components are now given by (using Equation (A6), (A13), and (A14))

(A18)

with U given by Equation (A14).

Comparing these expressions with those for the isotropic case (Equations (A1)), it follows that the anisotropic constitutive relation reduces to Glen’s law for α = β = 1 and γ = 2, with the exponent m = (n + 1)/2.

Application to plane flow

In the plane-flow approximation

and

must vanish, which implies τ′ _xy = τ′ _yz , = 0. Furthermore, because

incompressibility requires that

. From Equations (A18i) and (A18iii), we find

(A19)

while, from Equation (A18ii)

(A20)

Because, in the force-budget calculations, stresses as functions of strain-rates are required, the inverse formulation of Equations (A18) needs to be derived. Using Equations (A19) and (A20), it follows that

(A21)

Substituting these expressions in Equation (A14) yields

and hence

(A22)

with

and

The above expressions can be readily incorporated into the force-balance equations. Although not discussed here, it is obvious how this is done by checking the equations given in Van der Veen (1989). In the calculations, the parameter γ, describing enhanced vertical shearing, is taken to be proportional to the enhancement factors calculated by Dahl-Jensen (1985).

Applicability to glacier ice

To be theoretically correct, the tensors used in the above derivation should be the second-symmetric Piola-Kirchhoff stress tensor, and the Green-Saint-Venant strain-rate tensor. These tensors are defined with respect to a reference configuration in the material, rather than with respect to a fixed Cartesian coordinate system as used here.

For example, the Green-Saint-Venant strain-rate gives the rate of deformation from a reference configuration, whereas the classical strain-rate, as used in glaciology, uses spatial velocity gradients (cf. Truesdell and Toupin, 1960). The two strain-rate tensors are equal only for strains close to the reference configuration; that is, for short times during which the deformation and rotation remain infinitesimal. The Piola-Kirchhoff stress tensors are used in the material description of problems in which the boundary of a body is deformed in a prescribed way, or is subject to assigned forces (cf. Truesdell and Toupin, 1960).

For polar ice sheets, a reasonable choice for the reference configuration (X _r, Y _r Z _r) is one that is linked to the preferred orientation of the c-axes. For instance, the Z _r-axis could be taken along' the direction of preferred orientation, while the X _r- and Y _r-axes are chosen as such to form a right-handed, orthogonal coordinate system. The fixed coordinate system is defined as usual, with the z-axis vertical. The Green-Saint-Venant strain-rate is then defined with respect to the (X _r, Y _r Z _r) system, while the classical strain-rate refers to the fixed (x,y,z) system. As long as both coordinate systems remain fixed with respect to each other, the two strain-rates are identical. This will be the case if the c-axis fabric retains a vertical orientation.

It appears that this may be a valid first model. At a number of sites, vertical clustering of c-axes has been observed (e.g. Dye 3, Greenland (Herron and others, 1985); Camp Century, Greenland (Herron and Langway, 1982); Byrd Station, Antarctica (Gow and Williamson, 1976)). As shown by the analysis of Alley (1988), the pattern of flow at these sites (parallel flow at Byrd Station and Dye 3; divergent flow at Dome C and Camp Century) causes the c-axes to rotate toward the vertical axis.

This result may not be applicable to other parts of ice sheets. Therefore, a complete theory describing deformation of anisotropic glacier ice should consist of a constitutive relation linking the commonly used Cauchy stress to the classical strain-rate, and an equation describing the change with time of the director of anistropy, n. Such theories have been developed by Green (1964a, b), and in a series of papers by Ericksen (see Truesdell and Noll (1965, p. 523) for a discussion of Ericksen’s theory, and for references to the original papers). In Green’s theory, n is determined each time by the rotation from the reference configuration as described by the axis of transverse isotropy. Ericksen, on the other hand, derived a differential equation from which n can be calculated at any time.

None of the available theories has been adapted to the flow of glacier ice (and doing so obviously falls outside the range of competence of the present authors). A reasonable first step is therefore to investigate theories for anisotropy that can explain the observed pattern of flow near Dye 3, for the current configuration. Subject to this restriction, it appears that the anisotropy law discussed above is the best available at this time.

Appendix B

Normal Stress Effects In The Flow Of Glaciers

Flow law I modified second-order fluid

The first relation proposed by Man and Sun (1987) reads:

(B1)

where A_ij ⁽¹⁾ and A_ij ⁽²⁾ are the first and second Rivlin-Ericksen tensors, respectively. These are defined as

(B2)

and

(B3)

where summation over repeated indices is implied and u_k represents the velocity component in the k-direction. The function U in Equation (B1) is the second invariant of the strain-rate tensor

(B4)

If steady-state creep can be assumed, the time-derivative in Equation (B3) can be omitted. Furthermore, in the plane-flow approximation, the only non-zero strain-rates are

and

Equation (B3) then reduce to (denoting the x- and z-components of velocity by u and w, respectively)

(B5)

The stress-deviator components are then:

(B6)

with

(B7)

Inspection of these expressions shows that, for arbitrary values of α₁ and α₂, the three normal deviatoric stresses do not sum to zero. This is not an artefact of the plane-flow assumption, but is a feature of the general constitutive relation (BÍ). After some tedious arithmetic, the sum of the normal deviatoric stresses yields

(B8)

By definition, this sum must be zero, and the two parameters α₁ and α₂ are related as

(B9)

To choose values for the remaining three flow parameters (μ, α₁ and m), compare the first part of flow law I,

with Glen’s law

To allow for comparison between these two constitutive relations, we use

(B10)

and

(B11)

Having no better alternative, the remaining flow parameter is chosen to be proportional to μ. That is

(B12)

and hence, from Equation (B9)

(B13)

in which A ₁ is a constant of initially, unknown value.

Flow law It power-law fluid of grade 2

The second constitutive relation discussed by Man and Sun (1987) reads:

(B14)

Using the expressions given above, the non-zero stress-deviator components are:

(B15)

Again, the two parameters α₁ and α₂ are related through Equation (B9). The other flow parameters are chosen as above (Equation (B10), (B11), and (B12)).

Appendix C

Are Calculated Patterns Of Basal Drag Real Or Numerical Artefacts?

Stability of the calculation scheme

A worrying issue with the solution scheme developed to calculate stresses and strain-rates throughout a section of glacier using measured surface data, is that it may not be stable and that calculated horizontal variations in basal drag may be the result of numerical instabilities, rather than be real.

Numerical instabilities are a well-known phenomenon to any numerical modeller, and arise because of the discrete nature of numerical models. For example, in a time-evolving ice-sheet model, all relevant quantities (thickness, velocity, etc.) are calculated at grid points covering the area of interest. Time marching is achieved by solving the prognostic continuity equation for the change in ice thickness using a discrete time step. The size of this time step cannot be chosen independently from the spacing of the numerical grid. The Courant-Friedrichs-Lewy (CFL) condition provides a criterion for keeping the time integration stable (that is, inhibit the development of unwanted, numerically induced variations that grow explosively with time). In short, this condition states that the closer spaced the numerical grid the smaller the time step must be (cf. Smith, 1985).

The force-budget-solution technique is comparable to a time-evolving model. Only, instead of marching forward with time, calculations progress downward into the ice. Thus, it may be expected that a criterion similar to the CFL condition exists for keeping the downward calculation stable. That is, the number of depth layers (or the thickness of each layer) must be chosen in accordance with the horizontal grid spacing. However, due to the complicated nature of the equations involved (Van der Veen and Whillans, 1989a), finding a stability criterion by analytical means is no trivial matter. Therefore, results of various calculations, using different horizontal spacings, are compared to determine whether the outcome (calculated pattern of basal drag) is real or a numerical artefact.

To this purpose, velocities and elevations along the Dye 3 strain network are gridded to a regular grid with a spacing of 0.455 km, using cubic-spline interpolation (the survey poles of the strain grid are spaced about 1.7 km apart). Values thus obtained are used as input in four calculations, differing only in the horizontal grid spacing used (0.455, 0.91, 1.82, and 2.73 km, respectively), omitting in-between values for the less densely spaced numerical grids.

Effect of horizontal grid spacing on calculated patterns of basal drag

The driving stress, shown in Figure 10, is directly linked to gradients in the surface-elevation profile. Because increasing the horizontal grid spacing effectively smooths this profile, the curves in this figure exhibit progressively less variation as the grid-point distance is increased.

Because plane flow is considered here, side shear is neglected and the only resistive forces are due to differential pushes and pulls, and drag at the glacier base. The latter is the result of the force-budget calculation, being the residual obtained from consideration of forces. At the surface, the pushing or pulling resistive force is proportional to the second derivative of the surface velocity (shown in Figure 11). Again, increasing the horizontal grid spacing results in smoother curves. Small-scale features disappear, and the amplitude of the remaining variations becomes smaller, the larger the grid spacing. This implies that for large grid-point distances, calculated basal drag will be closer to the driving stress than when a small grid spacing is used (because of the smaller variations in longitudinal stretching).

Fig. 10. Driving stress. Values are obtained by gridding measured surface elevations to a regular grid with a spacing of 0.455 km, using cubic-spline interpolation. These values are assigned to four numerical grids with increasing spacing (0.455, 0.91, 1.82, and 2.73 km, respectively, from top to bottom), by omitting in-between values. These resulting surface elevations are used to calculate the driving stress. Tick marks on the horizontal bars under each curve indicate positions of the numerical grid points.

Fig. 11. Second derivative of the surface velocity (

). Values are obtained similarly as for Figure 10.

Fig. 12. Calculated patterns of basal drag from four calculations differing only in the horizontal grid spacing used. Value of the rate factor used is A₀ = 2 × 10⁵ kPa⁻³ year⁻¹, n = 3, and an Arrhenius-type temperature-dependence with the temperature as shown in Figure 3. The conventional generalized flow law is used and the enhancement factor equals 1 at all depths.

Calculated basal drag patterns (Fig. 12) confirm this. Horizontal variations become much smaller when the grid spacing is increased. For the largest spacing (2.73 km), calculated basal drag is virtually the same as the driving stress. For the second largest grid spacing (1.82 km), some differences between basal drag and driving stress occur (notably, the second maximum in basal drag at about 15 km is shifted somewhat to the right compared to the corresponding peak in driving stress) but, on the whole, both stresses follow a very similar pattern.

Further reduction of the grid-point distance (to 0.91 and 0.455 km) results in important differences between calculated basal drag and driving stress. Although the second curve shown in Figure 12 has a large-scale pattern corresponding to that of the two curves below it (with the exception of the minimum around 5 km), small-scale variations, not obvious in the corresponding curves in Figures 10 and 11, become important. For the smallest grid spacing used, the short wavelength features completely dominate the pattern of basal drag, and the resulting curve bears little resemblance to the other curves shown in Figure 12.

At first glance, these results may suggest that the numerical solution scheme becomes unstable for small horizontal grid sizes. However, this is not the case, as can be checked by increasing the number of depth layers in the calculation. Increasing this number from 21 (used in the calculations on which Figure 12 is based) to 101, and choosing the smallest horizontal spacing (0.455 km), yields virtually the same results as those represented by the upper curve in Figure 12 (differences between values of basal drag are only a few kPa, at most). To draw the analogy to time-marching models again, suppose a time-integration over 1000 years is carried out. If the results are the same when using a large time step of, say, 100 years, as they are for a much smaller time step, for instance 10 years, the results are considered to be real, and not caused by numerical instabilities.

The pattern of basal drag depicted in the upper curve of Figure 12 is not realistic. In places, drag at the base of the glacier is negative, tha.t is, directed down-flow. This artificial result is largely attributable to the short wavelength features. By applying a Gaussian filter to the calculated pattern of basal drag, the results become much more realistic (Fig. 13). The broad-scale pattern corresponds reasonably well with that calculated using the larger horizontal grid spacing (1.82 km). Differences are still apparent but these can be traced back to differences in the input data (Figs 10 and 11).

Origin of small-scale variations in basal drag

Because the short wavelength features dominating the pattern of basal drag calculated using the smallest grid spacing are not caused by numerical instabilities, there must be a physical explanation for their existence.

It is well established that surface features on an ice sheet reflect basal conditions. For example, a site of high basal drag may result in a locally steep surface, or a velocity variation at the surface, or a combination of the two. However, the ice acts as a filter so that the amplitude of horizontal variations decreases going upward. In other words, basal variations show up at the surface, but with a muted amplitude. The extent of this glacial filtering (decrease in amplitude) is strongly dependent on the wavelength of the features under consideration. The filter function (ratio of surface to basal amplitude) rapidly goes to zero as the wavelength decreases to less than about one ice thickness (e.g. Hutter and others, 1981; Whillans and Johnsen, 1983). Thus, a basal feature with a horizontal scale of, say, half the ice thickness, will be difficult to detect at the surface if nresent at all.

Fig. 13. Effect of Gaussian smoothing on the pattern of basal drag calculated using the smallest horizontal grid spacing (0.455 km). Panel a shows the unsmoothed pattern, as in the upper panel of Figure 12; panel b that was obtained after Gaussian smoothing with a band width of twice the grid spacing, and panel c after smoothing with a band width of four times the grid spacing. Panel d shows the calculated basal drag from the calculation using a horizontal grid spacing of 1.82 km, as in the third panel of Figure 12.

The force-budget calculation uses measured surface velocities and slopes to calculate basal drag. Hence, it is an inverse calculation, in which observed effects are used to determine the causes responsible for these effects. The filter works now in the opposite way, to amplify surface effects as the calculation proceeds downward and short wavelength features are most amplified. If the horizontal grid spacing is sufficiently small, this selective amplification may result in short wavelength variations dominating the calculated pattern of basal drag.

It should be remarked here that the solution technique does not involve simplifying assumptions or approximations, other than that of plane flow (which in this case, is not severe restriction). Given that the numerical scheme appears to be stable, this implies that calculated basal drag is as realistic as the input data are.

The short-scale variations in Figure 12 are therefore attributed to deficiencies in the input data for the calculations. The deficiencies arise from the interpolation scheme used to grid the measured data and from normal measurement uncertainties. It is reassuring that smoothing the results of a too-detailed calculation (Fig. 13) results in about the same basal drag variation obtained from a more appropriate calculation.