Skip to main content Accessibility help


  • Access



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Mass Balance Along Two Transects of the West Side of the Greenland Ice Sheet
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Mass Balance Along Two Transects of the West Side of the Greenland Ice Sheet
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Mass Balance Along Two Transects of the West Side of the Greenland Ice Sheet
        Available formats
Export citation


The mass balance is computed along the Ohio State University (OSU) transect near the Arctic Circle and along the Expédition Glaciologique Internationale au Groenland (EGIG) line. Measured surface velocities are compared with velocities calculated from up-glacial accumulation rate, flow-line spreading, ice thickness, and the depth variation in horizontal velocity. The depth variation in velocity is calculated using the constitutive relation for ice, calculated temperatures within the glacier, computed shear and longitudinal stresses, and allowance for impurity content and ice-crystal orientation. The resulting mass balance is +0.6 ± 0.14 m a−1 for the OSU transect and 0 ± 0.07 m a−1 along the EGIG line. The errors arise mainly from uncertainties in measured accumulation rate and flow-line spreading, and perhaps in flow-enhancement factors due to ice anisotropy or impurities. The results for the EGIG line differ from prior estimates mainly because earlier works placed greater emphasis on short-term accumulation rates.


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 Bindschadler (1984). The OSU (Ohio State University) transect west from the south dome has not been previously analyzed, except for a brief summary by Kostecka and Whillans (1986).

Data along the EGIG transect were obtained by surface surveying for motion (Hofmann, 1975). Accumulation rates are from pole heights measured over a 9 year period (Mälzer and Seckel, 1976) and calculated by Radok and others (1982), and from interpretation of firn strata (Mock and Weeks, 1965; Reeh and others, 1978). Ice thicknesses are by airborne radar (personal communication from S. Overgaard), near-surface temperatures from Mock and Weeks (1965), and elevation contours from Bindschadler (1984).

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


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 Zwally 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 ui 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 (Whillans and Cassidy, 1983; van 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 uz 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 (Doake and Wolff, 1985; Wolff 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 [ux (bed) = 0]

for the horizontal velocity, ux (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 (Budd, 1968), along the Byrd Station Strain Network (Whillans and Johnsen, 1983), and near Dye 3, Greenland (Whillans 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 −10 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−1 at temperatures less than −10°C and 139 kJ mole−1 at higher temperatures (Paterson, 1981) are used. The base value of A −10 enters only in calculating stretching stress, and the value of 5.2 × 10−16 s−2 kPa−3, the softness parameter at −10°C (Paterson, 1981, p. 39), is used. In the equation for shape factor, constant proportional uncertainties in A −10 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 Whillans (1981)). No major temperature changes seem to have occurred since the Holocene warming (Dansgaard 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 (Yen, 1981, table 3):

for θ in K. C y is volumetric thermal capacity (Bolzan, 1984, p. 6):

and Q is internal heat production (Paterson, 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−2 using data from southern Greenland (Sass 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 (Raymond, 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 Bolzan (1985), Budd (1969, p. 48–53), Hammer and others (1978, p. 22), Whillans (1979), and Whillans and Cassidy (1983). Lliboutry (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 (Whillans, 1979), and Camp Century, Greenland (Budd 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 Doake 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 −10 (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 Raymond (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−1 for the OSU transect and 0.02 m a−1 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.


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 (Zwally and others, 1983; Bindschadler, 1984). These directions agree within about 5° of the velocity vectors measured on the OSU transect (Drew 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 (Hofmann, 1975) are corrected as suggested by Robin (1983, fig. 2.17b) by subtracting a supposed systematic error of about 3.5 m a−1 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 Zwally 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 (Weidick, 1975) and satellite radar (Bindschadler, 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 Mock (1967) found, some factor in addition to elevation is important. There is no relation with local slope (Whillans 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 (Radok and others. 1982). and the circles are long-term averages from deep-core stratigraphy (Reeh 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 (Whillans 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. Reeh 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−2 a−1, or 15% (1 standard deviation) (Whillans, 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−1 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 Reeh and others (1978) and of Reeh 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−2 to 0.0505 W m−2 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. Budd and others, 1976; Whillans, 1979; McInnes 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 Shoji and Langway (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 (Gundestrup 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 Budd and others (1976). McInnes and Radok (1985). and Whillans (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 Nye (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 (Whillans 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. Shoji and Langway (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 −10, 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 −10 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 Gundestrup 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.

Fisher 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. Wolff 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.


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−1). Thus, the mass balance of the OSU transect is +0.08 − (5 ± 20)% × 0.4 = +0.06 ± 0.08 m a −1 or between −0.02 and +0.14 m a −1. Along the EGIG transect it is 0.02 − (5 ± 17)% × 0.4 = 0 ± 0.07 m a−1.


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. Mälzer and Seckel (1976) compared surface elevations in 1959 and 1968, and found a mean thickening of 0.06—0.07 m a −1 which is consistent with our results. Just to the south, Bindschadler (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 −4 a−1 for the model with no ice thickening. The measured total strain-rate there is 1.23 × 10−4 a−1 (Reeh 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.

Grigoryan and others (1985) proposed that the ice sheet is currently thickening at 0.2 m a−1. 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 (Whillans, 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 (Gundestrup and Hansen, 1984) became available, it could be argued that uncertainty in the enhancement factors of Shoji and Langway (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 (Zwally, 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.


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.


Bindschadler, R.A. 1984 Jakobshavns glacier drainage basin: a balance assessment. Journal of Geophysical Research 89(C2), 206672.
Bolzan, J.F. 1984 Ice dynamics at Dome C, East Antarctica. Ohio Slate University. Institute of Polar Studies. Report 85.
Bolzan, J.F. 1985 Ice flow at the Dome C ice divide based on a deep temperature profile. Journal of Geophysical Research, 90(D5), 811124.
Budd, W. 1968 The longitudinal velocity profile of large ice masses. International Association of Scientific Hydrology Publication 79 (General Assembly of Bern 1967 — Snow and Ice), 5877.
Budd, W.F. 1969 The dynamics of ice masses. ANARE Scientific Reports, Ser. A(IV), Glaciology. (Publication 108.)
Budd, W.F. Young, N.W. 1983 Application of modelling techniques to measured profiles of temperatures and isotopes. In Robin de, G.Q. ed. The climatic record in polar ice sheets. Cambridge, etc., Cambridge University Press, 15077.
Budd, W.F. Young, N.W. Austin, C.R. 1976 Measured and computed termperature distributions in the Law Dome ice cap, Antarctica. Journal of Glaciology 16(74), 99110.
Dansgaard, W. Clausen, H.B. Gundestrup, N. Johnsen, S.J. Rygner, C. 1985 Dating and interpretation of two deep Greenland ice cores. In Langway, C.C.jr, Oeschger, H., Dansgaard, W. eds. Greenland ice core: geophysics, geochemistry, and the environment. Washington, DC, American Geophysical Union, 7176. (Geophysical Monograph 33.)
Doake, C.S.M. Wolff, E.W. 1985 Flow law for ice in polar ice sheets. Nature, 314(6008), 25557.
Drew, A.R. Whillans, I.M. 1984 Measurement of surface deformation of the Greenland ice sheet by satellite tracking. Annals of Glaciology, 5, 5155.
Fisher, D.A. Koerner, R.M. 1986 On the special rheological properties of ancient microparticle–laden Northern Hemisphere ice as derived from bore–hole and core measurements. Journal of Glaciology, 32(112) 50110.
Grigoryan, S.S. Buyanov, S.A. Krass, M.S. Shumskiy, P.A. 1985 The mathematical model of ice sheets and the calculation of the evolution of the Greenland ice sheet. Journal of Glaciology, 31(109), 28192.
Gundestrup, N.S. Hansen, B.L. 1984 Bore–hole survey at Dye 3, south Greenland. Journal of Glaciology 30(106), 28288.
Hammer, C.V. Clausen, H.B. Dansgaard, W. Gundestrup, N. Johnsen, S.J. Reeh, N. 1978 Dating of Greenland ice cores by flow models, isotopes, volcanic debris, and continental dust. Journal of Glaciology 20(82), 326.
Heeswijk van, M. 1984 Meteorite concentration by ice flow. Ohio State University. Institute of Polar Studies. Report 83.
Hofmann, W. 1975 Die Internationale Glaziologische Grönlandexpedition (EGIG). 2. Die geodätische Lagemessung. — Eisbewegung 1959–1967 in den EGIG–Profilen. Zeitschrift für Gletscherkunde und Glazialgeologie, 10, 1974 21724.
Kostecka, J.M. Unpublished. A numerical ice flow model applied to mass balance calculations in western Greenland. (M.Sc. thesis, Ohio State University, 1985)
Kostecka, J.M. Whillans, I.M. 1986 Mass balance along two transects of the Greenland ice sheet. (Abstract.) Annals of Glaciology, 8, 198.
Lliboutry, L. 1981 A critical review of analytical approximate solutions for steady state velocities and temperatures in cold ice–sheets. Zeitschrift für Gletscherkunde und Glazialgeologie, 15(2), 1979 13548.
McInnes, B. Radok, U. 1985 A steady–state prediction of Dye 3 core features. In Langway, C.C.jr, Oeschger, H., Dansgaard, W. eds. Greenland ice core: geophysics, geochemistry, and the environment. Washington, DC, American Geophysical Union, 11117. (Geophysical Monograph 33.)
Mälzer, H. Seckel, H. 1976 Internationale Glaziologische Grönlandexpedition (EGIG). Zeitschrift für Gletscherkunde und Glazialgeologie, 11(2), 1975 24552.
Mock, S.J. 1967 Calculated patterns of accumulation on the Greenland ice sheet. Journal of Glaciology, 6(48), 795803.
Mock, S.J. Weeks, W.F. 1965 The distribution of ten–meter snow temperatures on the Greenland ice sheet. CRREL Research Report 170.
Nye, J.F. 1959 The motion of ice sheets and glaciers. Journal of Glaciology, 3(26), 493507.
Paterson, W.S.B. 1981 The physics of glaciers. Second edition. Oxford, etc., Pergamon Press.(Pergamon International Library.)
Radok, U. Barry, R.G. Jenssen, D. Keen, R.A. Kiladis, G.N. McInnes, B. 1982 Climatic and physical characteristics of the Greenland ice sheet. Parts I and II. Boulder, CO, University of Colorado. Cooperative Institute for Research in Environmental Sciences.
Raymond, C.F. 1983 Deformation in the vicinity of ice divides. Journal of Glaciology, 29(103), 35773.
Reeh, N. Clausen, H.B. Gundestrup, N. Johnsen, S.J. Stauffer, B. 1977 δ (18O) and accumulation rate distribution in the Dye 3 area, south Greenland. International Association of Hydrological Sciences Publication 118 (General Assembly of Grenoble 1975 — Isotopes and Impurities in Snow and Ice), 17781.
Reeh, N. Clausen, H.B. Dansgaard, W. Gundestrup, N. Hammer, C.U. Johnsen, S.J. 1978 Secular trends of accumulation rates at three Greenland stations. Journal of Glaciology, 20(82), 2730.
Robin de, G.Q. 1983 The Greenland ice sheet. In Robin de, G.Q. ed. The climatic record in polar ice sheets. Cambridge, etc., Cambridge University Press, 3842.
Sass, J.H. Neilsen, B.L. Wollenberg, H.A. Munroe, R.J. 1972 Heat flow and surface radioactivity at two sites in south Greenland. Journal of Geophysical Research, 77(32), 643544.
Shoji, H. Langway, C.C. 1984 Flow behavior of basal ice as related to modeling considerations. Annals of Glaciology, 5, 14148.
Weidick, A. 1975 A review of Quaternary investigations in Greenland. Ohio State University. Institute of Polar Studies. Report 55.
Whillans, I.M. 1976 Radio–echo layers and the recent stability of the West Antarctic ice sheet. Nature, 264(5582), 15255.
Whillans, I.M. 1977 The equation of continuity and its application to the ice sheet near “Byrd” Station, Antarctica. Journal of Glaciology, 18(80), 35971.
Whillans, I.M. 1978 Surface mass–balance variability near “Byrd” Station, Antarctica, and its importance to ice core stratigraphy. Journal of Glaciology, 20(83), 30110.
Whillans, I.M. 1979 Ice flow along the Byrd Station strain network, Antarctica. Journal of Glaciology, 24(90), 1528.
Whillans, I.M. 1981 Reaction of the accumulation zone portions of glaciers to climatic change. Journal of Geophysical Research, 86 (C5), 427482.
Whillans, I.M. Cassidy, W.A. 1983 Catch a falling star: meteorites and old ice. Science, 222(4619), 5557.
Whillans, I.M. Johnsen, S.J. 1983 Longitudinal variations in glacial flow: theory and test using data from the Byrd Station strain network, Antarctica. Journal of Glaciology, 29(101), 7897.
Whillans, I.M. Jezek, K.C. Drew, A.R. Gundestrup, N. 1984 Ice flow leading to the deep core hole at Dye 3, Greenland. Annals of Glaciology, 5, 18590.
Whillans, I.M. 1987 Glaciologic transect in southern Greenland: 1980 and 1981 GISP–I data. >Columbus, OH, Ohio State University. Byrd Polar Research Center. (Report 2.)
Wolff, E.W. Doake, C.S.M. 1986 Implications of the form of the flow law for vertical velocity and age–depth profiles in polar ice. Journal of Glaciology, 32(112), 36670.
Yen, Y. 1981 Review of thermal properties of snow, ice and sea ice. CRREL Report 8110.
Zwally, H.J. 1986 Ice–sheet thickening observed by satellite altimetry. (Abstract.) Annals of Glaciology, 8, 200.
Zwally, H.J. Bindschadler, R.A. Brenner, A.C. Martin, T.V. Thomas, R.H. 1983 Surface elevation contours of Greenland and Antarctic ice sheets. Journal of Geophysical Research, 88(C3), 158996.

Appendix Calculation Scheme

1. Use smoothed surface and bed topographies, flow line spreading R(x), thickness H(x), accumulation rate b(x), [m(x) = 0], secular rate of thickness change specified (initially zero), and calculate mean horizontal “continuity” velocity, u x (x), using

where x represents distance along flow line and t, time.

2. At the ice divide, assume vertical velocity varies linearly with depth:

and calculate steady-state temperature profile, θ(z), at ice divide

with K = K(θ) and C v = C v(θ).

This is used as an initial condition for the temperature calculation of the next step.

3. Make first estimate of depth variation in horizontal velocity:

The value of p in this first estimate is not critical; we use p = 2. Elevation of top surface is represented by z s.

4. Use velocities calculated from

and vertical velocity from incompressibility:

with boundary conditions at the top surface (which is more efficient than at the bed). Compule steady state temperatures with depth and distance along flow line:

5. Solve for horizontal-stretching stress,


6. Calculate an improved shape factor:

Go to step 4.

7. Solution accepted on convergence of