Introduction
A re-analysis of bore-hole tilt data from four sites has recently been carried out (Doake and Wolff, 1985). The data suggested that a linear flow law may be as appropriate as the traditional power law for interpreting the flow of ice in polar areas. Clearly, this would have important implications for many features of ice sheets. As part of an attempt to determine past precipitation rates, Reference Paterson and WaddingtonPaterson and Waddington (1984) have discussed vertical velocity profiles and age–depth profiles, using a power law with n = 3. Here, we consider the effect on similar calculations of using a linear-flow law (including temperature in our modelling). The two problems to be investigated are (i) profiles of velocity with depth at the Devon Island Ice Cap, and (ii) the age–depth profile at Camp Century, Greenland.
Theory
Much of the relevant theory has been discussed by Reference RaymondRaymond (1983) and Reference Paterson and WaddingtonPaterson and Waddington (1984).
The flow law for ice is generally expressed as
where έ is effective strain-rate, τ effective stress, n the flow-law exponent, and A is a factor that is dependent both on temperature and on properties of the ice such as fabric and impurity content. A is often expressed as
where R is the gas constant, T is temperature, and Q, the activation energy, is often taken to be 60 kJ mol–1 (Reference PatersonPaterson, 1981), a value used throughout this paper.
We now adopt a coordinate system in which distance x and velocity u are horizontal, z and w are vertical (positive upwards) and z = 0 at the base of the ice. Then, from laminar-flow theory,
where έ xz is shear strain-rate, h is ice thickness, α is the surface slope, ρ is the density of ice, and both z and h are expressed in ice-equivalent depths throughout. Integrating from the base of the glacier, the horizontal velocity profile with depth is given by:
where u s is the measured velocity at the surface. Equation (4) assumes that A 0 is constant with depth.
Using the condition of incompressibilily, allowing either u s or α and h to be slowly varying functions of x, and neglecting a small term in (Reference RaymondRaymond, 1983; Reference BolzanBolzan, 1984)
Again, integrating from the bed, the vertical velocity profile with depth is given by
where w s is the surface vertical velocity, which, assuming that the ice thickness is not changing with time and that there is no bottom melting, can be taken as the ice-equivalent annual accumulation. In practice, f(z) and g(z) may have to be determined by summation rather than integration.
If we assume an isothermal ice sheet. Equations (4) and (6) can be determined by integration to give (Reference RaymondRaymond, 1983)
These formulae will not apply near an ice divide. Reference RaymondRaymond (1983) used Equations (7) and (8) as boundary conditions far from a divide. He then used finite-element modelling to derive u(z) and w(z) between the boundary and the divide. The results for n = 3 were shown in figure 3 of Reference RaymondRaymond (1983). For n = 1, the profiles are not sensitive to distance from the divide and conform to the laminar flow Equations (7) and (8) for isothermal, and Equations (4) and (6) for non-isothermal.
Finally, to derive an age-depth profile,
where t is age in years and λ is the annual layer thickness, both at height z above the bed. If the ice sheet is in steady state from the surface down to height z, this becomes the formula given by Reference Paterson and WaddingtonPaterson and Waddington (1984), which should read
in our notation.
We note that strictly there is an inconsistency in using Equation (4) or (7) for the horizontal velocity profile and Equation (6) or (8) for the vertical profile. The horizontal velocity profile is for the case of laminar flow in which there is no dependence of any quantity with x and would give a zero vertical velocity everywhere. A first-order correction to allow either u s or α and h to vary slowly with x has been made by Reference RaymondRaymond (1983) and Bolzan (1984) in their respective models to give a vertical velocity profile with depth and we have followed their example. A more rigorous treatment of the subject has been given by Reference Morland and JohnsonMorland and Johnson (1980). We indicate in the following sections the likely effects of the approximate nature of our analyses.
Vertical Velocity Profile at Devon Island
Both horizontal and vertical velocity profiles have been measured in a bore hole on Devon Island in Arctic Canada. The bore hole is only 2–3 ice thicknesses from the ice divide, so laminar-flow models cannot be used in the case n = 3. Reference Paterson and WaddingtonPaterson and Waddington (1984) adapted Reference RaymondRaymond’s (1983) finite-element model, developed for isothermal ice in plane strain, to the conditions at Devon Island, to compute profiles for n = 3. They made an allowance for temperature and had to make both a width factor W(x) (introduced to allow transverse deformation) and the flow-law factor A 0 adjustable parameters in order to fit the data. The resultant profiles showed an excellent fit to the measured vertical velocity profile but a poor fit to the horizontal velocity profile.
Since, for n =1, profile shape is not sensitive to distance from the divide, we can use Equations (4) and (6) to obtain profiles for n = 1, with temperatures from Reference RobinRobin (1983). The results of this, together with those for n = 3 given by Reference Paterson and WaddingtonPaterson and Waddington’s (1984) model, are shown in Figure 1. As we would have predicted from our earlier work (Doake and Wolff, 1985) on the stress–strain -rate relationship, n = 1 gives a much better fit to the horizontal velocity profile than does n = 3. Our computed vertical velocity profile for n = 1 is a poor fit to the data plotted in figure 4 of Reference Paterson and WaddingtonPaterson and Waddington (1984). Their n = 3 model fits these data well. However, the data plotted appear to be derived from far more scattered original data given by Reference PatersonPaterson (1976). In the original data, the error bar on each vertical velocity point was given as 0.032 m a–1. A large part of this error bar arises from the measurement of the surface vertical velocity, which acts as a systematic error on the individual measured velocities. Therefore, both the n = 1 and n= 3 models are well within the error bars.
We conclude that the data for the two profiles together are at least as well fitted by n = 1 as by n = 3, provided temperature is properly accounted for. Results from Reference RaymondRaymond (1983) suggest that the concavity of the depth profile of w might depend as much on longitudinal variations in the strain-rate (έ x ) as on the ice rheology. This indicates that a more rigorous approach may be necessary. However, the variation in longitudinal strain-rate near the bore hole on Devon Island is small (Reference PatersonPaterson, 1976) and likely to have an insignificant effect.
The present-day velocity profiles, unlike age–depth profiles, are not affected by past changes in the ice sheet such as changes in slope, ice thickness, or accumulation. They will be affected by present–day spatial variations in the structure of the ice. For instance, if A 0 varies with depth, then Equation (4) must be modified. It is quite possible that the A 0 of Wisconsin ice is different from that of Holocene ice. However, at Devon Island the Wisconsin corresponds to only the bottom few metres of ice, so even radically different values of A 0 in this ice can do little to alter the overall shape of the velocity profiles.
Age-Depth Profile at Camp Century
The Camp Century bore hole was drilled to bedrock. The site is far from the main ice divide, although some authors have speculated about the possible effects of a local divide (Reference Hammer, Hammer, Dansgaard, Gundestrup, Johnsen and ReehHammer and others, 1978). A time-scale has been established by various methods (Reference Hammer, Hammer, Dansgaard, Gundestrup, Johnsen and ReehHammer and others, 1978). Reference Paterson and WaddingtonPaterson and Waddington (1984) used Equations (8) and (10) to derive a flow-law time-scale for the past 8500 years. They used n = 3 and a number of assumptions, one of which was that the ice was isothermal (it is in fact ~ 11 K warmer at the base than the surface). Their time-scale was indistinguishable from that of Reference Hammer, Hammer, Dansgaard, Gundestrup, Johnsen and ReehHammer and others (1978), and this allowed them to make an assessment of past precipitation rates at this site.
We have calculated two further flow-law time-scales, one for n = 1 and one for n = 3. In both we use Q = 60 kJ mol–1, and then use summation in Equations (4), (6) and (10) to compute profiles. The ice-equivalent depth, h, is 1367 m. Following Reference Paterson and WaddingtonPaterson and Waddington (1984), w s is taken as 0.403 m of ice per year. The temperature profile is from Reference RobinRobin (1983).
Figure 2 shows that our n = 1 profile is indistinguishably close to Reference Hammer, Hammer, Dansgaard, Gundestrup, Johnsen and ReehHammer and others’ (1978) time-scale back to 8500 years B.P. Reference Hammer, Hammer, Dansgaard, Gundestrup, Johnsen and ReehHammer and others (1978) stated that 10 000 years occurs 236 m above the base, i.e. at 1131m. Our n = 1 time-scale would predict this age to be at 1136 m.
The n = 3 non-isothermal time-scale is well removed from Reference Hammer, Hammer, Dansgaard, Gundestrup, Johnsen and ReehHammer and others’ (1978) time-scale, and predicts 10 000 years at 1176 m. At 1000 m depth, the discrepancy is 40 m or 700 years (~ 10%), well outside the likely error (estimated at better than 3% over 10 000 years) in Reference Hammer, Hammer, Dansgaard, Gundestrup, Johnsen and ReehHammer and others’ (1978) time-scale. The agreement between the non-isothermal flow-law time-scale for n = 1 and the time-scale derived by Reference Hammer, Hammer, Dansgaard, Gundestrup, Johnsen and ReehHammer and others (1978) by independent methods suggests that the depth profile of w is relatively unaffected by variations in the longitudinal strain-rate (έ x ) with x. However, more data on έx would be needed to verify this conclusion. Beyond 10 000 B.P., both of our time-scales, in common with all other flow-law time-scales, predict ages that are far too young. This is assumed to be because there may have been a different flow regime and accumulation in Wisconsin times.
Several assumptions have been made concerning the structure and history of the ice sheet (as opposed to its physics) in our analysis. First, in moving from Equation (3) to the equation for the horizontal velocity, we assumed that A 0, α and Q were constant with depth. However, different depths may correspond to surface slopes averaged over a different horizontal scale, and A 0 and Q may depend on properties of the ice (chemistry, fabric, etc.) which change with depth.
Data analysis (using bore-hole tilting data with n = 3) by Reference PatersonPaterson (1983) has previously suggested that the Wisconsin-age ice at Camp Century may be considerably (at least ten times) harder than Holocene-age ice (i.e. A 0 is lower in Wisconsin ice). If this is so, then A 0 should be included in top and bottom of Equation (4). This would change the shape of the profiles of u and w. It would also alter the predicted age of younger ice as well as that of the Wisconsin ice itself. For n = 1, no significant change in A 0 with depth was required to fit bore-hole tilting data (Doake and Wolff, 1985). We have not therefore attempted to include A 0 changes in our calculations.
The second major assumption is that involved in going from Equation (9) to Equation (10). Here it is assumed that the vertical velocity is equal to λ, the annual layer thickness, down to height z for which the age, t, is calculated. This implies that the accumulation rate has been constant over this time period (and is the same at the sites further up-stream where the snow fell). It also implies that all the ice down to height z has been subjected to the same thinning processes. This in turn requires that α and A 0 have not varied significantly and have in the past been the same at all depths. These conditions are, of course, not expected to hold for Wisconsin ice and may account for the failure of any model to predict correct ages for Wisconsin ice. However, if this assumption fails only below height z, it does not affect the age at heights above z. Past changes in the ice sheet in the Wisconsin do not therefore invalidate our model for Holocene ages, provided that the Holocene ice is in steady state (and that A 0 is currently constant with depth, or its variation is explicitly modelled). A more general approach would be to integrate along a particle path, thus following variation in distance as well as depth.
A final, and potentially serious, problem at Camp Century is that it is rather close to a local ice divide. If this has an important effect on the flow, then this, and any other, analysis of behaviour at Camp Century is invalid, and any empirical conclusions will be of purely local interest.
Although trends in these variables can invalidate the model, “local” changes will not, but will merely cause fluctuations in layer thickness about the values that would be predicted. Reference Hammer, Hammer, Dansgaard, Gundestrup, Johnsen and ReehHammer and others (1978) measured layer thicknesses at 22 points in the top 1100 m at Camp Century. If all of our other assumptions hold, but short-term variations in accumulation rate, a t, are allowed, then
Using w s/w calculated from the models, and λ as measured by Reference Hammer, Hammer, Dansgaard, Gundestrup, Johnsen and ReehHammer and others (1978), we can then calculated at. Using this method for n = 3 isothermal, Reference Paterson and WaddingtonPaterson and Waddington (1984) suggested that short-term variations in accumulation rate of ±20% have occurred over the last 8500 years at Camp Century. Our n = 1 temperature-dependent model naturally gives a similar result.
It is also possible that a particular layer of ice may have (and always have had) a different A 0 from the remaining ice, due to a different fabric, chemistry, etc. Layers with a higher A 0 would have thinned faster than the rest, and this could be falsely interpreted as evidence that a t was lower than the mean. To investigate the likely scale of this effect, we can hold a t constant, but allow A 0 to vary about a mean value (while remaining on average constant). Vertical strain-rate is defined by
where σz’ is the vertical stress-deviator component.
Then, for a layer,
where λ is the layer thickness at time t, and a 0 was the original thickness of the layer ai its accumulation date.
Now, consider two ice sheets, identical except that the values of A 0 in a thin layer at depth z are (and always have been) different. The integral on the right of Equation (13) is therefore the same in the two ice sheets. Then:
If we now set A 01 equal to the mean A 0 in the rest of the ice sheet, and λ1 equal to the λ predicted from our model, then using the measured λ values we can calculate the true A 0 in each layer. We have used this method to calculate A 0 rel = A 0/A0 average for Reference Hammer, Hammer, Dansgaard, Gundestrup, Johnsen and ReehHammer and others’ (1978) data, assuming a 0 constant. The variation of A0 rel with depth is shown in Figure 3. Except at the very top, A 0 need only vary by ±15% to produce the observed λ pattern. It has been suggested that the A 0 value of Wisconsin ice at Camp Century could be at least ten times lower than that of Holocene ice (for n = 3) (Reference PatersonPaterson, 1983) or, alternatively, that the A 0 value of Wisconsin ice at other sites may be 3–4 times greater than that of Holocene ice (Reference ReehReeh, 1985). Although it is not yet clear which, if either, of these views is correct, it seems reasonable to suggest that smaller changes in A 0 (of up to 15%) could have occurred within the Holocene, due to changes in ice chemistry or fabric.
Therefore, although fluctuations in accumulation rate undoubtedly occur, it is not possible to be conclusive about past accumulation rates from layer thicknesses unless we know how A 0 varies with depth. This requires a very accurate knowledge of changes in strain-rate with depth. It also suggests that more research is needed to determine which ice parameters affect the value of A 0.
Discussion and Conclusions
We have discussed two problems in terms of the flow law of polar ice. For Devon Island, we have previously shown indirectly (Doake and Wolff, 1985) that n = 1 fits the horizontal velocity profile, while n = 3 does not (Reference Paterson and WaddingtonPaterson and Waddington, 1984). For the vertical velocity profile, both flow laws give predicted results well within the error bars of the measurement.
More dramatic results are found in assessing the age–depth profile at Camp Century. It had previously been shown that an n = 3 isothermal model fitted the accepted time-scale very well. However, we have shown that a more realistic non-isothermal (Q = 60 kJ mol–1) model for n = 3 fits very poorly. A non-isothermal model for n = 1 gives a near-perfect fit. The non-isothermal n = 3 model can be forced closer to the data by using a lower w s (by up to 12%) but the overall shape is still a poor fit.
We therefore conclude that, at Camp Century, as with the tilt data (Doake and Wolff, 1985), the age–depth data suggest that n = 1 may give a more appropriate description of behaviour. Models with n = 3 can be made to fit only if Q = 0 (which gives results equivalent to the isothermal case), if the accumulation decreases dramatically and monotonically with depth, or if the flow-law constant A 0 decreases significantly in deeper ice. If the local ice divide is important, then the n = 3 models at least are not valid. We also note that, whatever flow law is used, any variations of A 0 with depth could be as important as changes in accumulation rate in determining annual layer thicknesses. Therefore, attempts to derive past accumulation rates from changes in layer thickness must at present be treated cautiously.
It appears that scenarios that fit the Camp Century bore-hole tilt data (linear-flow law, or n = 3 with decreasing A 0 with depth) also fit the age–depth profile to 10 000 years. Therefore, putting the measured tilt data into our equations would also lead to a good fit to the age–depth profile to 10 000 years. The bore-hole tilt data and associated measurements at Camp Century are not considered very accurate, and cannot by themselves define a flow law. However, the age–depth profile supports the tilt data in the sense that the same models fit both.
In the two problems we have studied here, we have shown that a linear law gives as good or better a fit than is obtained with n = 3. On the other hand, it is perhaps surprising that the assumption of perfect plasticity (Reference ReehReeh, 1982) gives a realistic reconstruction of the shape of the central Greenland ice sheet, though elevations are underestimated in places. Ice-sheet profiles are determined largely by deformation near the base. This result indicates that, with a suitable choice of basal parameters, a range of flow exponents can reproduce the essential features of the ice sheet. However, in order to model in detail the flow of ice sheets, past and present, accurate flow-law parameters are still required. These are likely to be obtained most easily from bore-hole data to as great a depth as possible. Many further situations might be considered. It is noteworthy that in one of these problems we found that n = 1 (Q = 60 kJ mol–1) and n = 3 isothermal gave almost identical profiles. In other words, temperature variations at this site were giving the same effect as an extra two powers of the flow-law exponent. As many ice-sheet features have, for mathematical convenience, been modelled isothermally with n = 3, we may expect similar data to emerge in some cases from non-isothermal models with n = 1. The Camp Century profile in particular shows that it is not satisfactory to assume that an ice sheet is isothermal; the results this gives can often be far removed from those obtained when temperature is included in the model.
Addendum
What is the value of n, the flow-law exponent?
The current paper does not attempt to answer this question, but rather to examine the effect on real data of using two different values. The data here are not sufficient to define a flow law, but they do add weight to the idea that the a priori assumption of n = 3 may sometimes lead to a poor fit to data obtained in the field.
Reference Jezek, Jezek, Alley and ThomasJezek and others (1985) have suggested that a good fit to n = 3 is obtained using data from bottom crevasses on the Ross Ice Shelf. A number of factors, however, make these data unsuitable for the determination of the flow law.
-
1. The state of stress in the ice shelf was derived from a theory relating height of crevasse to total ice thickness (Reference WeertmanWeertman, 1980). However, it is uncertain whether or not the crevasse heights as measured by radio echo-sounding are in equilibrium with the in-situ stress field.
-
2. The theory of crevasse formation on the underside of ice shelves has not been tested by other data sets. The assumption that the theory is perfect in the form in which it was originally presented (Reference WeertmanWeertman, 1980) weakens the case for placing trust in results derived from applying the theory to a situation where there are additional uncertainties.
-
3. The variation of temperature within the ice shelf has not been allowed for. The difference in temperature between the ice front and the grounding line is sufficient to cause the strain-rates to vary by a factor of three.
-
4. Strain-rates were not measured at the sites where crevasse data have been used to infer the stress. Instead, strain-rates have been interpolated from a regular grid with an average spacing of about 50 km between nearest points.
-
5. There are discrepancies between the stress fields in the ice shelf calculated in one case from bottom crevasse data (Reference Jezek, Jezek, Alley and ThomasJezek and others, 1985) and in another case from using measured strain-rates in an assumed power flow law (Reference Thomas and MacAyealThomas and MacAyeal, 1982). In both cases, the same form of flow law was used (or derived) and the same strain-rate data set used. Thus the difference between the two sets of results (sometimes exceeding 105 Pa) should be an indication of the possible error in the effective stress parameter. Errors of this magnitude render unsuitable any data set for the purpose of deriving a realistic flow law.
If the factors above were properly accounted for, they would tend to push the highly scattered data towards a lower value of n.
High-quality bore-hole tilt data are potentially available from the bore hole at Dye 3. However, all publications on this to date have been based on the differences between two surveys. The first of these was carried out at the time of drilling, and used an inclinometer in which “readings … tend to be too low, in particular at inclinations approaching the upper limit of the inclinometer (6° )” (Reference Gundestrup and HansenGundestrup and Hansen, 1984). Unfortunately, the inclinometer reached 5° at two-thirds depth, and was within a fraction of 6° for the bottom 500 m. It is therefore hard to see how any conclusion can be reached about the flow at Dye 3 from these data. Other authors (Reference Dahl-JensenDahl-Jensen, 1985; Reference Shoji and LangwayShoji and Langway, 1985) have used the data without comment. Further inclinometer surveys at this site will be eagerly awaited.
It might be concluded that the best that can be said of attempts to define the flow law is that the data used so far have all been of poor quality, and the analyses carried out with them may therefore be dismissed. This perhaps indicates merely that the field does not provide the ideal conditions required either for models to be deduced or for them to be applied. However, it also suggests that it is desirable that authors should consider both of the postulated flow laws, as we have done here, in studies of polar ice sheets.