Introduction

A challenge in calculating mass balance from measured surface velocities and up-glacial accumulation is estimating by how much the surface velocity exceeds the mean velocity through the thickness. This difference can be obtained from data on the progressive tilting of a bore hole, but there are few bore holes. Alternatively, as done here, the correction can be determined from calculation of the major deviatoric stresses and the constitutive relation. The result is then compared with measured tilting. If the test is successful, the model can be applied elsewhere where there are no bore holes for the computation of mass balance.

The calculation, however, requires the constitutive relation for ice (which describes how easily the ice deforms under a given stress), which is not well known. Simple calculations using stress balance and taking account of thermal effects and crystal-orientation fabric produce a surface velocity that is accurate to a factor of two, at best. And so simple calculations of internal shear are inadequate in calculating mass imbalances, which are likely to be less than 50% of the current discharge. Much of the uncertainty arises within the deepest layers where stresses and temperatures are both large, and small uncertainties in quantities in this zone of great shear affect the velocities through the entire thickness.

Much of this difficulty can be avoided by using the constitutive relation to calculate not the total velocity but the difference between surface and mean velocity. This difference is typically only 10–30% of the surface velocity and the contributions of the problematic deep layers are much smaller. Thus, if only the difference between surface and mean velocity is sought, the constitutive relation can be used with much greater confidence.

This approach is taken here and the velocity variation with depth is described by a shape function that is normalized to the mean velocity. In this way, the determination of mean velocity is separated from the problem of velocity variation with depth. The mean velocity is computed using the equation of continuity with specified secular thickening or thinning along the flow line. The shape function is computed using this mean velocity together with calculations of temperature and longitudinal stress, and perhaps including fabric enhancement. The result is a predicted value of surface velocity for the specified mass balance (rate of thickening or thinning). Comparison with measured velocities indicates whether the mass balance was correctly specified.

The separation of the determination of mean velocities from the determination of the velocity—shape function means that most of the uncertainty in the value of the softness parameter in the constitutive relation is unimportant. This scheme requires only relative depth variations in softness parameter due to variations in temperature or structure.

The two studied transects are shown in Figure 1. The EGIG (Expédition Glaciologique Internationale au Groenland) transect has been studied in this regard before, most recently by Reference BindschadlerBindschadler (1984). The OSU (Ohio State University) transect west from the south dome has not been previously analyzed, except for a brief summary by Reference Kostecka and WhillansKostecka and Whillans (1986).

Data along the EGIG transect were obtained by surface surveying for motion (Reference HofmannHofmann, 1975). Accumulation rates are from pole heights measured over a 9 year period (Reference Mälzer and SeckelMälzer and Seckel, 1976) and calculated by Reference Radok, Barry, Jenssen, Keen, Kiladis and McInnesRadok and others (1982), and from interpretation of firn strata (Reference Mock and WeeksMock and Weeks, 1965; Reference Reeh, Clausen, Dansgaard, Gundestrup, Hammer and JohnsenReeh and others, 1978). Ice thicknesses are by airborne radar (personal communication from S. Overgaard), near-surface temperatures from Reference Mock and WeeksMock and Weeks (1965), and elevation contours from Reference BindschadlerBindschadler (1984).

Data for the OSU transect were obtained by Doppler satellite tracking and short-arc orbit calculations for precise positions (Reference Drew and WhillansDrew and Whillans, 1984), and detection of the nuclear bomb levels in firn (Reference WhillansWhillans and others, 1987) and earlier stratigraphic interpretation (Reference Mock and WeeksMock and Weeks, 1965). Ice thicknesses were obtained by airborne radar (Reference WhillansWhillans and others, 1987; personal communication from S. Overgaard). Temperatures were measured at 17 m depth (Reference WhillansWhillans and others, 1987) and flow-line spreading was obtained both from the measured velocity field and from radio-echo mapping by satellite (Reference Zwally, Bindschadler, Brenner, Martin and ThomasZwally and others, 1983).

Theory

Fig.1. Studied transects in southern Greenland. The EGIG transect runs west from Crête and the OSU transect west from near Dye 3. Base map with elevation contours (m) is simplified from Reference Zwally, Bindschadler, Brenner, Martin and ThomasZwally and others (1983).

The first part of the calculation provides a velocity based on the equation of continuity:

in which H represents ice thickness, b and m surface net accumulation and basal melting rates in terms of ice of mean glacial density, t time, and x and y horizontal coordinates. The mean velocity component in the i-direction (i = x, y) is:

where z represents the vertical coordinate, positive upward, and u_i the horizontal component of velocity in the i-direction.

This continuity relation is simplified by defining the x-axis to follow the flow line. The y-component of velocity, u _y is then zero along the x-axis. Transverse spreading, d u _y/dy, is described by a divergence parameter, R, which, because flow is usually perpendicular to elevation contours, is also the radius of curvature of those contours (Reference Whillans and JohnsenWhillans and Cassidy, 1983; Reference Heeswijk vanvan Heeswijk, 1984):

Adjustments due to the flow line being curved are negligible. Applying these relations to the equation of continuity, provides

along the flow line (y = 0).

This equation is solved numerically for the mean horizontal velocity, u _x (x) using data for H(x), b (x), R (x), and prescribed values for . Basal melting is set to zero because later temperature calculations show that the basal ice is frozen to the bed along both flow lines ( Fig.3). The position of the ice divide, where u _x =0 is specified.

The depth variation in horizontal velocity, , is obtained from the constitutive relation (third power):

Here u_z represents vertical velocity, A is the flow-law “softness” parameter, and horizontal shear stress is represented by τxy. This is the conventional flow law. There are reservations about the flow law (Reference Doake and WolffDoake and Wolff, 1985; Reference Wolff and DoakeWolff and Doake, 1986), but this form is in general usage and other forms were not tried. The effect of different exponents is discussed in the section on sensitivity. Horizontal gradients in vertical velocity, , are ignored and we find a posteriori that to be a valid approximation. The effective shear stress τ_e, is given by:

where are normal deviatoric stresses. In this work, τ _yz is zero because the x-axis follows the flow line. For the flow lines studied, the transverse stresses, τ _xy and are also very small and, although is included in the continuity calculation, both it and the terms and are associated with very small stresses that are neglected here. Thus, because normal deviatoric stresses sum to zero,

and the constitutive relation used is

which is integrated from the bed [u_x (bed) = 0]

for the horizontal velocity, u_x (x, z).

It is the ratio of horizontal velocity to the mean velocity that is needed here and we call it the “shape factor”, ψ. That is:

Of special note is that because the stresses and softness parameter, A, appear in similar fashions in denominator and numerator, depth-independent uncertainties in them cancel and so are not important. Only relative depth variations in stresses and softness parameter affect the shape factor. The two relevant stresses and the softness parameter are discussed, in turn, next.

The horizontal shear stress, τ _xz , balances the driving stress if gradients in longitudinal stress can be neglected Thus,

is used for horizontal shear stress, in which ρ represents mean ice density, g the magnitude of the acceleration due to gravity, and z _s the elevation of the ice-sheet surface. Measured gradients in longitudinal strain-rate on Law Dome, Antarctica (Reference BuddBudd, 1968), along the Byrd Station Strain Network (Reference Whillans and JohnsenWhillans and Johnsen, 1983), and near Dye 3, Greenland (Reference Whillans, Jezek, Drew and GundestrupWhillans and others, 1984) vary importantly only when scales of less than 15 km are considered. The need to consider longitudinal gradients in stress is therefore avoided by binomially smoothing the surface and bed topographies at a scale of 7 km (one standard deviation). The horizontal shear stress described by the above equation is therefore appropriate for the smoothed flow.

Although gradients in normal stress, , may on average be small, the normal stress, , itself in general is not. Moreover, at shallow depths and near the ice divide where the shear stress, τ _xz , is small, this normal stress is proportionately more important. It is calculated within an iterative loop (Appendix, steps 4, 5, and 6) from the velocities and shear stress using the constitutive relation for stretching:

which is a cubic equation in with only one real root. Figure 2 shows results for two positions in Greenland and it is evident that normal stresses are very important at shallow depths and near the ice divide.

Fig.2. Stress distribution through thickness. τ _xz is horizontal shear stress. is longitudinal normal deviator. and τ _e is the effective shear stress calculated from these.

The stretching stress softens the ice to horizontal shear because of its role in the effective shear stress, and we slightly underestimate its importance. This is because we use smoothed surface and bed topographies which leads to much reduced variation in stretching stress. Neglect of this effect leads to calculated horizontal shearing at shallow depths that is too small. This issue is discussed further under “sensitivity”.

The flow-law parameter depends on ice temperature and mechanical properties according to:

in which A ₋₁₀ is the flow-law parameter at −10°C (263K.), depending on crystal anisotropy and impurity content. The gas constant is represented by R _g and θ(x, z) is ice temperature in K. For the activation energy, Q, values of 60 kJ mole⁻¹ at temperatures less than −10°C and 139 kJ mole⁻¹ at higher temperatures (Reference PatersonPaterson, 1981) are used. The base value of A ₋₁₀ enters only in calculating stretching stress, and the value of 5.2 × 10⁻¹⁶ s⁻² kPa⁻³, the softness parameter at −10°C (Reference PatersonPaterson, 1981, p. 39), is used. In the equation for shape factor, constant proportional uncertainties in A ₋₁₀ cancel.

The calculation of softness parameter, A, and horizontal shear requires temperature and that is a function of both depth and distance along the flow line. Although the continuity velocity calculated earlier allows for non-steady conditions , non-steady temperature effects are not considered. As it turns out, the results are not critically sensitive to temperature and so this is not an important limitation. Also, the temperature wave associated with the Holocene warming 10 000 years ago penetrated into the bed 5000 years ago for the OSU transect, and along the EGIG line adjustment should be largely complete more recently (using the theory in Reference WhillansWhillans (1981)). No major temperature changes seem to have occurred since the Holocene warming (Reference Dansgaard, Clausen, Gundestrup, Johnsen, Rygner, Langway, Oeschger and DansgaardDansgaard and others, 1985). Thus, neglect of secular temperature change is also not an important limitation.

Heat transfer is described by:

in which horizontal conduction is neglected and K represents ice conductivity (Reference YenYen, 1981, table 3):

for θ in K. C _y is volumetric thermal capacity (Reference BolzanBolzan, 1984, p. 6):

and Q is internal heat production (Reference PatersonPaterson, 1981, p. 200):

The small work done by stretching and lateral divergence is neglected. Boundary conditions are measured surface temperature as a function of position, and basal heat flux which is taken to be 0.0431 W m⁻² using data from southern Greenland (Reference Sass, Neilsen, Wollenberg and MunroeSass and others, 1972). The temperature calculation is done by finite differences using the initial conditions at the ice divide described in the Appendix. A separate solution is made for each iteration on the velocity field. Results are shown in Figure 3.

Fig.3. Calculated isotherms for 5°C intervals.

At the ice divide, horizontal velocities are zero and the scheme for calculating vertical velocities in step 4 does not work. Instead, a simple estimate of the depth variation in vertical velocity is made (linear with depth). This is not accurate (Reference RaymondRaymond, 1983) and so results near the ice divide must be treated with caution.

A first estimate for the velocity—depth shape function is used following Reference BolzanBolzan (1985), Reference BuddBudd (1969, p. 48–53), Reference Hammer, Clausen, Dansgaard, Gundestrup, Johnsen and ReehHammer and others (1978, p. 22), Reference WhillansWhillans (1979), and Reference Whillans and CassidyWhillans and Cassidy (1983). Reference LliboutryLliboutry (1981) provided some theoretical justification for these parameterizations. The versions vary slightly but they are all similar to:

in which z _s is the elevation of the top surface and ξ is the proportion of surface velocity that is due to internal shear (as opposed to basal sliding). The multiplying fraction in the definition of ψ ensures that the thickness mean of ψ equals 1. In Greenland there is no basal sliding and ξ = 1. The exponent (p + 1) is obtained by empirical fits to measured or calculated velocity profiles, and Figure 4 compares versions of Equation (2) with calculated shape functions.

Fig.4. Calculated shape function, ψ versus height above bed for two sites along the EGlG transect (solid lines), together with examples of ψ according to Equation (2) (dashed lines).

The form of the velocity—depth profile is interesting in itself. If the ice were isotropic, isothermal, and deformed only by horizontal shear, then p in Equation (2) would be the exponent, 3, in the flow law. In practice, higher temperature at depth and crystal anisotropy make p larger, and the effect of normal stresses makes p smaller. Fits to the velocity profile at Byrd Station, Antarctica (Reference WhillansWhillans, 1979), and Camp Century, Greenland (Reference Budd, Young and Robin deBudd and Young, 1983, fig. 5.42), show p = 1 and 1.5, respectively. This suggests that normal stresses are very important at those sites, or alternatively, as Reference Doake and WolffDoake and Wolff (1985) argued, that the exponent in the flow is small, but we use the conventional flow law. At great distances from ice divides our work ( Fig.4) shows that the flow is more nearly constant with depth (indicated by larger values of p) because longitudinal stresses are proportionately less important.

The next sections discuss the results more fully but the main features are shown in Figures 2—6, which are for the case of constant softness parameter, A ₋₁₀ (no allowance for anisotropy or impurity content but temperature effects included). The depth variation in stresses at two sites along the EGIG line is shown in Figure 2 and the important contribution of longitudinal stress to the effective shear stress is evident. The profiles are similar to those calculated by Reference RaymondRaymond (1983) near a hypothetical ice divide. The associated velocity-profile shape functions, ψ, are shown in Figure 4. The surface value of ψ, [ψ (z _s)], is the scaled difference between surface velocity and the mean through the thickness; it is used in the mass-balance calculation and its variation along the flow lines is shown in Figure 5. Beyond 50 km from the ice divide, surface velocity is 10—20% larger than the mean velocity through the thickness, and the difference is greater near the ice divide because of the reduced proportional contribution of horizontal shear stress to the effective shear stress. These values of ψ(z _s) are used to calculate steady-state surface velocities and they are found to be slightly larger than those measured ( Fig.6). The lower panels in Figure 6 show that the velocities match for the cases of secular thickening of 0.08 m a⁻¹ for the OSU transect and 0.02 m a⁻¹ for the EGIG line. Uncertainties are, however, important and these are the subject of a sensitivity study.

Fig.5. Surface values of ψ for both flow lines

Fig.6. The dots represent measured velocities and are the same for the upper and lower figures. The lines are calculated velocities.

Sensitivity

Parameters that are comparatively well constrained are thickness and flow-line direction. Thicknesses are from airborne radio echo-sounding. Flow lines used here are drawn perpendicularly to surface-elevation contours obtained by analysis of satellite-borne radar altimetry (Reference Zwally, Bindschadler, Brenner, Martin and ThomasZwally and others, 1983; Reference BindschadlerBindschadler, 1984). These directions agree within about 5° of the velocity vectors measured on the OSU transect (Reference Drew and WhillansDrew and Whillans, 1984). The precise locations of flow lines are not important because thickness measurements show no strong across-flow gradients, and so the calculations here follow radio-echo flight lines which approximate the flow directions determined by the methods above.

Calculated velocities are compared with measured values. The measured values from along the EGIG transect (Reference HofmannHofmann, 1975) are corrected as suggested by Reference Robin de and Robin deRobin (1983, fig. 2.17b) by subtracting a supposed systematic error of about 3.5 m a⁻¹ in the north—south direction. Possible errors in the correction scheme are not critical. Along the OSU transect the velocities are well determined by repeated Doppler-satellite tracking with short-arc orbits calculated from stationary tracking stations on the western coast of Greenland. In comparison with other aspects, the measured velocities on the OSU transect can be considered almost error-free.

Lateral spreading or convergence can be measured with strain networks or calculated from the radius of curvature of elevation contours, R. Along the OSU transect these methods can be compared. Radii of curvature on the map of Reference Zwally, Bindschadler, Brenner, Martin and ThomasZwally and others (1983) are +200 km and −480 km (480 km is large and indicates nearly parallel flow) for the central and western clusters (shown in Figure 1). The spreading indicated by measured velocities shows some scatter, because locally the flow turns a little in association with basal obstructions; however, mean values for R from measured velocities are 175 km and near infinity (parallel flow) for the central and western clusters, respectively. Thus, there is close agreement and the maximum difference between the “satellite-radar” map and the measured spreading leads to an uncertainty of only 4% in velocity calculated using continuity.

Along the EGlG transect the elevation contours obtained by surface surveying (Reference WeidickWeidick, 1975) and satellite radar (Reference BindschadlerBindschadler, 1984) show that the flow is nearly parallel and the radius of curvature is set to near infinity in the calculations. Errors in this are difficult to assess because large-scale lateral strain-rates are not available for that area. The uncertainty in calculated surface velocity due to lateral spreading is assumed to be similar to that along the OSU transect (4%).

Accumulation rates are very critical and affect calculated surface velocities proportionately. There is a large variability and a geographic pattern in accumulation rate is not evident and interpolation between measurements is uncertain. Measured accumulation rates along the OSU transect are shown in Figure 7 and, as Reference MockMock (1967) found, some factor in addition to elevation is important. There is no relation with local slope (Reference WhillansWhillans and others, 1987). The problem along the EGIG transect is, as discussed below, that different techniques produce different values ( Fig.7), but the range of values is less than along the OSU transect. The figures include straight-line representations of the accumulation rate—elevation relations used in the calculations.

Fig.7. Accumulation rale against elevation. The open triangles are from stake measurements for a 9 year interval (Reference Radok, Barry, Jenssen, Keen, Kiladis and McInnesRadok and others. 1982). and the circles are long-term averages from deep-core stratigraphy (Reference Reeh, Clausen, Dansgaard, Gundestrup, Hammer and JohnsenReeh and others. 1978). The diamonds are from pit studies of unstated time interval (Mock and Weeks. 1965). The solid triangles are accumulation rates for the interval 1965—80 from detection of nuclear-bomb fall-out levels (Reference WhillansWhillans and others. 1987). The center lines are the favored relations used in this study and the other lines represent uncertainties.

The three straight-line relations are drawn to include most of the range of values, placing greater reliance on longer-term average values. (We know of no physical model that would support a statistical fit.) The central lines are used in the main calculations. The uncertainties represented by the other lines are 12% for the OSU transect and 9% for the EGIG transect and, as Table I shows, they are the major source of error in the mass-balance calculation.

There is a systematic bias between measurement techniques for accumulation-rate data along the EGIG transect. Stake measurements show accumulation rates about 18% less than does pit work. Core studies lie in the center of the range. The core results are for centuries-long time averages, and the pit and stake results may indicate short-term secular variation. Reference Reeh, Clausen, Dansgaard, Gundestrup, Hammer and JohnsenReeh and others (1978) indicated that 30 year means can deviate by about 10% from the long-term mean in this area. Shorter-term means would deviate even more. Thus, secular variation can account for the differences and the secular variance is important to mass-balance studies.

For comparison, in central West Antarctica annual variability is 20 kg m⁻² a⁻¹, or 15% (1 standard deviation) (Reference WhillansWhillans, 1978). The pit and stake measurements in Greenland span 4 years (estimated, because Mock and Weeks did not describe the data source) and 9 years, respectively, and the standard deviation of variability for these time averages would be 7.5% and 5% or about 0.025 m a⁻¹ by analogy with central Antarctica. This is about the same as the observed differences between accumulation rates for different time spans. Thus, the bias between measurement techniques is reasonably accounted for by inter-annual variability.

The relation with elevation along the OSU transect is weak ( Fig.7). Large variation is also observed in the data of Reference Reeh, Clausen, Dansgaard, Gundestrup, Hammer and JohnsenReeh and others (1978) and of Reference Reeh, Clausen, Gundestrup, Johnsen and StaufferReeh and others (1977) from their work near Dye 3 and the local ice divide where two 5 year mean accumulation rates differ by 40%. The greater variation may be because the OSU transect is farther south and closer to the open ocean and so maritime influences are correspondingly greater. There may be additional explanations, but the scatter in Figure 7 means that well-constrained accumulation rates cannot be obtained, and that is a severe limitation to the mass-balance calculations.

Uncertainties in the temperature model are not critical. This is because, unlike a full stress and strain-rate calculation, the temperatures are used only to compute the velocity-shape factor, ψ, whose thickness mean equals 1 by definition. Calculated surface velocities are the product of mean velocities, obtained from continuity, and the surface value of ψ, and so it is ultimately only ψ(z _s that is needed for the mass balance.

The most poorly constrained parameter in the temperature calculation ( Fig.3) is the basal heat flux. As an example, changing it from the favored 0.0431 W m⁻² to 0.0505 W m⁻² leads to an average increase in basal temperature of 2° C and an average decrease of surface velocity of 2%. The decrease occurs because in each calculation the mass balance is constrained, and warmer basal ice means that a greater proportion of the shearing occurs at depth and that correspondingly less occurs at shallow depths. The velocity profile is then more plug-like and the surface velocity nearer in value to the mean (as represented by ψ = 1). Reasonable variations in other aspects of the temperature calculation show even smaller effects on the surface velocity.

Fitting the simple parameterization (Equation (2)) to the calculated velocity profile (as shown in Figure 8) suggests p ≃ 8 for no sliding (ξ= 1). This value of p is much larger than that inferred from bore-hole tilting at Camp Century and Byrd Station, but does agree with that inferred from the bore-hole tilting at Dye 3 ( Fig.8). We have no explanation for this major contrast between Dye 3 and the other core holes.

Vertical velocities are calculated from incompressibility and horizontal strain-rates (step 4 in the Appendix), and these can be compared with suggested approximations for variation of the vertical velocity with depth (e.g. Reference Budd, Young and AustinBudd and others, 1976; Reference WhillansWhillans, 1979; Reference McInnes, Radok, Langway, Oeschger and DansgaardMcInnes and Radok, 1985). Our calculated vertical-velocity profiles ( Fig.9) are not, in general, similar to the simple parameterizations. This is because the vertical velocity is strongly affected by basal slope and can even change sign, as at 80 km on the OSU transect, because of reverse slopes. This finding is very relevant to calculations of the depth—age relationship for ice cores. They are obtained from integrals of vertical-velocity profiles like those shown in Figure 7, and there are important differences between the parameterizations and the velocities calculated here for different sites.

Fig.8. Enhancement factors versus relative height for Dye 3 and Camp Century from Reference Shoji and LangwayShoji and Reference Dansgaard, Clausen, Gundestrup, Johnsen, Rygner, Langway, Oeschger and DansgaardLangway (1984). and the effect of their inclusion on the velocity-profile shape factor at 20 km from the ice divide on the EClG transect. Pluses are unenhanced shape factors for 40 km on OSU transect, corresponding to the distance of Dye 3 from ice divide, on the other side. Dashed line is the measured profile at Dye 3 (Reference Gundestrup and HansenGundestrup and Hansen, 1984).

Fig.9. Vertical velocity through the thickness at three positions along the OSU transect. Both velocity and height are normalized to surface values. The solid lines are calculated and the dotted lines are simple parameterizations that have been suggested by Reference Budd, Young and AustinBudd and others (1976). Reference McInnes, Radok, Langway, Oeschger and DansgaardMcInnes and Radok (1985). and Reference WhillansWhillans (1979) for the case of no basal sliding (ξ = 1).

It is important to include longitudinal stress in the calculation because it softens the ice to horizontal shear through the effective shear stress. Figure 2 shows example depth distributions of horizontal shear stress, longitudinal stress, and effective shear stress. As found by Reference NyeNye (1959) near the ice divide, the resulting effective shear stress is, to a first order, constant with depth. This is in great contrast to lamellar (or laminar) flow models in which the viscosity varies as the square of depth due to stress softening. Here, the viscosity varies with depth mainly due to other effects (temperature and ice-flow enhancement). Neglect of longitudinal stress causes calculated surface velocities to be 30% smaller near the ice divide and 10% smaller away from the divide.

Although horizontally averaged longitudinal stresses are included, the model does not include the effects of local variations. Stresses are calculated for a binomially smoothed glacial geometry in order to eliminate the complication of longitudinal-stress gradients in calculating horizontal shear stress. However, the longitudinal-stress variations mean that locally the longitudinal stress is larger or smaller than the smoothed, mean value. Along the Byrd Station Strain Network in Antarctica, local longitudinal stretching ranges from about twice the regional value to a small fraction of the mean (Reference Whillans and CassidyWhillans and Johnsen, 1983). For longitudinal strain-rate varying between zero and twice the value used in the smoothed model, the effective stress at the glacial surface would vary between its value for no contribution by longitudinal stresses and (2)^½ times larger (because the stretching enters as a square in step 6 of the Appendix). The correction to surface velocity as calculated here is −10% where stretching is zero and +(2)^½ 10% where stretching is large. Thus, on average, the smoothed model should be corrected for the effect of longitudinal stresses by ½((2)^½ − 1) 10% or 2%. In Table I, surface velocities are increased by this amount.

The final important parameter is the enhancement of ice shear due to ice anisotropy and impurity content. Reference Shoji and LangwayShoji and Reference Dansgaard, Clausen, Gundestrup, Johnsen, Rygner, Langway, Oeschger and DansgaardLangway (1984) have proposed regions of the Greenland ice sheet where strain-rate is as much as 12 times the value otherwise calculated. That is, the softness parameter, A ₋₁₀, should be as much as 12 times larger than used in the standard calculation at the indicated depths. As with the temperature calculations, the depth variation in A ₋₁₀ does not affect the mean velocity, only the shape factor and the ratio of surface to mean velocity changes. Inclusion of enhancement factors decreases calculated surface velocity by 7%.

Enhancement factors are not included in the basic calculation because the accuracy of the enhancement factors is not known. They were determined from laboratory tests at unnaturally high stress, and their variation along the flow line is completely unknown. However, using the present model, the enhancement factors of Shoji and Langway do predict a velocity profile ( Fig.8) that roughly matches that measured at Dye 3 by Reference Gundestrup and HansenGundestrup and Hansen (1984). Thus, we believe them to be good estimates and, based on the fit to the Dye 3 profile, assign an uncertainty of 2% to calculated surface velocities due to uncertainties in enhancement factors.

Reference Fisher and KoernerFisher and Koerner (1986) argued for enhancement factors of only about 3. In that case, the adjustment for enhancement would be —2%.

Calculations are performed for a third power constitutive relation (n = 3, in the usual notation), but other suggested power relations are not expected to change the results very much. The power affects the exponent of the effective shear stress in Equation (1), being 2 in this case. Reference Wolff and DoakeWolff and Doake (1986) suggested that this exponent should be 0 (n = 1) and that would remove the effect of depth-varying effective shear stress. However, as is evident from Figure 2, near the ice divide the effective shear stress is nearly constant with depth and changing the value of n has the same effect as a change in softness parameter which, as noted above, is unimportant here. Farther from the ice divide, the effective shear stress shows more variation with depth ( Fig.2) such that a smaller exponent decreases the relative amount of shearing at depth and increases the difference between the surface and mean velocity, and leads to a more positive calculated mass balance. However, because of the important role of longitudinal stretching stress, especially at shallow depths, the longitudinal-velocity profile and our calculations are not greatly sensitive to the exponent in the flow law.

Conclusion

Table 1 summarizes the adjustments to the calculations displayed in the figures. The net adjustment is −5% and error limits are 20% for the OSU transect and 17% for the EGIG line. The effect on the calculated mass balance is obtained by multiplying by the mean accumulation rate (0.4 m a⁻¹). Thus, the mass balance of the OSU transect is +0.08 − (5 ± 20)% × 0.4 = +0.06 ± 0.08 m a ⁻¹ or between −0.02 and +0.14 m a ⁻¹. Along the EGIG transect it is 0.02 − (5 ± 17)% × 0.4 = 0 ± 0.07 m a⁻¹.

Discussion

The objective has been to calculate the mass balance and to study the limitations to its accuracy. Conventional methods are used and the commonly recognized effects are included. In particular, temperature effects and longitudinal stresses are included in calculating the difference between surface and mean velocities. Both studied transects have mass balances that are not significantly different from zero and the main uncertainty is due to scatter in measured values for the accumulation rate.

The result for the EGIG transect agrees with earlier work. Reference Mälzer and SeckelMälzer and Seckel (1976) compared surface elevations in 1959 and 1968, and found a mean thickening of 0.06—0.07 m a ⁻¹ which is consistent with our results. Just to the south, Reference BindschadlerBindschadler (1984) compared drainage through Jakobshavns Fjord with accumulation rate and found thickening but not significantly different from zero. The accumulation rates were determined by pit interpretation and, as shown in Figure 7, these values are larger for the EGIG line than others and that probably explains Bindschadler’s calculated thickening. At the ice divide, the calculated surface stretching is 1.3 × 10 ⁻⁴ a⁻¹ for the model with no ice thickening. The measured total strain-rate there is 1.23 × 10⁻⁴ a⁻¹ (Reference Reeh, Clausen, Dansgaard, Gundestrup, Hammer and JohnsenReeh and others, 1978), which is in good agreement considering the uncertainty in our calculation of strain-rate near an ice divide. Thus, all studies based on field data show central Greenland in near balance.

Reference Grigoryan, Buyanov, Krass and ShumskiyGrigoryan and others (1985) proposed that the ice sheet is currently thickening at 0.2 m a⁻¹. This is at variance with the present results and we call into question the importance of the simplifications in that mathematical model.

The main limitation to accuracy is the observed scatter in accumulation rate. Because of this, long-term values for accumulation rate cannot be estimated with precision. This problem contrasts with similar work near Byrd Station, Antarctica (Reference Drew and WhillansWhillans, 1977), where the accumulation rate was well determined. The problem is that accumulation rate in Greenland is more irregularly distributed in position or time, or both. That may be because of the more maritime climate of Greenland, and more precise mass-balance calculations by the present method require long-term averages of accumulation rate such as can be obtained from long ice cores.

Until the results of the Dye 3 core-hole tilting (Reference Dansgaard, Clausen, Gundestrup, Johnsen, Rygner, Langway, Oeschger and DansgaardGundestrup and Reference Gundestrup and HansenHansen, 1984) became available, it could be argued that uncertainty in the enhancement factors of Reference Shoji and LangwayShoji and Reference Dansgaard, Clausen, Gundestrup, Johnsen, Rygner, Langway, Oeschger and DansgaardLangway (1984) is a major limitation. However, because those enhancement factors together with allowance for longitudinal stress and temperature correctly predict the measured core-hole tilting ( Fig.8), more confidence is placed in those enhancement factors.

These problems highlight the advantages of repeat radar or laser altimetry by satellite (Reference ZwallyZwally, 1986). The issues raised here do not enter into the calculation of mass balance by the satellite techniques; they enter rather in the interpretative stage when the long-term significance and cause of any thickness change is discussed.

Finally, we remark on a part of the study that is not fully discussed in this paper. The calculations were extended to determine particle trajectories and isochrones (Kostecka, unpublished). The shapes of the isochrones were then compared with internal radar layers along the EGIG line. Within measurement and model uncertainties the two sets agree. However, because there is not a strong horizontal gradient in accumulation rate, neither radar layers nor isochrones show a distinctive tilt. Without such an accumulation-rate gradient, many non-steady models of ice-sheet flow yield similar isochrones and the techniques cannot be used to infer past mass balance, as Whillans (1976) was able to do in West Antarctica.

Acknowledgements

We thank R. Alley, J. Bolzan, and C.J. van der Veen for discussions and advice. Drafting is by R. Tope. This work was supported by U.S. National Science Foundation grant DPP 8008356. This is Byrd Polar Research Center contribution No. 601.