Skip to main content Accessibility help


  • Access
  • Cited by 114


      • 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.

        Force Budget: I. Theory and Numerical Methods
        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.

        Force Budget: I. Theory and Numerical Methods
        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.

        Force Budget: I. Theory and Numerical Methods
        Available formats
Export citation


A practical method is developed for calculating stresses and velocities at depth using field measurements of the geometry and surface velocity of glaciers. To do this, it is convenient to partition full stresses into lithostatic and resistive components. The horizontal gradient in vertically integrated lithostatic stress is the driving stress and it describes the horizontal action of gravity. The horizontal resistive stress gradients describe the reactions. Resistive stresses are simply related to deviatoric stresses and hence to strain-rates through a constitutive relation.

A numerical scheme can be used to calculate stresses and velocities from surface velocities and slope, and from ice thickness. There is no mathematical requirement that the variations in these quantities be small.


The flow of glaciers is driven by gravity and opposed by resistive forces such that there is zero net force acting on any section of a glacier (inertial effects being negligible). The gravitational driving force is well understood but a major problem in understanding glacial flow is the proper identification and recognition of the sites of action of the resistive forces. These can act at the bed (basal drag), at the sides of a glacier or flow band (side drag), and longitudinally as pushes or pulls. The calculation of these resistive forces is an important component of the study of the flow and stability of glaciers.

This has been a long-standing problem in glaciology. Nye (1957) calculated the shear stress at depth for the case of very simple bed configurations and plane flow. This approach was followed and improved upon by Robin (1967), Collins (1968), Nye (1969), Budd (1969, 1970a, b), and Hutter (1983), who incorporated the effects of gradients in longitudinal pushes and pulls. These studies show that basal resistive drag, τb, is linked to the driving stress, −⍴gH∂h/∂x with a correction for longitudinal gradients and a “variational” stress term, commonly denoted by T:

Here, x is directed along the flow and represents the vertical mean longitudinal deviatoric stress. To calculate it is often assumed that the longitudinal strain-rate is independent of depth (Nye, 1957; MacAyeal and others, 1987).

The T-term on the right-hand side of the above expression for basal drag is rather complex, involving the second derivative of the shear stress σxz,:

in which z represents the vertical coordinate. In most applications, the T-term is small and may be omitted from force-budget calculations (Hutter, 1983; Whillans and Johnsen, 1983; Kamb and Echelmeyer, 1986). The physical implications of omitting it are, however, not conceptually straightforward.

Further improvements in force-balance studies include side drag through a shape factor (e.g. Raymond, 1980), or through another gradient term (Raymond, 1980; Bindschadler and others, 1987; Frolich and others, 1987).

All previous investigations on the stress fields in glaciers rely on simplifying assumptions. For instance, the otherwise rigorous analysis of Hutter and others (1981) is limited to cases in which the bottom and top surfaces deviate only slightly from straight parallel lines. Observations of real glaciers indicate, however, that the basal relief is often very large (10% or more of the ice thickness, over a horizontal distance of few kilometers). Similarly, the theory of Kamb and Echelmeyer (1986) is based on a linearization in which longitudinal variations in the ice flow are treated as small perturbations to the overall average flow. Also, longitudinal stresses are supposed not to contribute significantly to the effective stress. Longitudinal stresses in ice sheets are, however, very large and may play important roles.

Another major problem in the practical application of previous theories is that these require basal conditions to be known. In general, however, conditions prevailing at the base of a glacier (basal drag and sliding velocities) are not known and, with most existing methods, it is necessary to vary these conditions until an agreement between, for example, calculated and measured surface velocities is obtained. This makes the method of solution somewhat arbitrary and time-consuming.

What is needed is a method for using measurements of the surface effects of flow variations to calculate stresses and strain-rates at depth. In particular, the need exists to use measured data to calculate resistive stresses opposing the ice flow, and to identify the sites of stress concentration.

Earlier attempts to develop such a method are those of Whillans and Johnsen (1983) and Whillans and Jezek (1987). In these studies, basal conditions are not prescribed, but resistive stresses are calculated from measured surface-velocity and slope data. These methods, however, cannot treat power-law creep and are restricted to plane flow.

In this paper, a generally applicable method for the interpretation of field data in terms of the stress field is developed. First, the equations describing balance of forces are derived using the partitioning of full stresses into lithostatic and resistive parts (as in Whillans (1987)). To relate these stresses to observable strain-rates, the constitutive relation is invoked and the balance equations are rewritten in terms of strain-rates. Approximate solutions can be readily obtained from these equations by making simplifying assumptions. Exact solutions require numerical techniques, which are explained in the second part of this contribution.

Deviatoric stresses are needed for the constitutive relation which links stresses to strain-rates. However, it is easier to understand force budget in terms of resistive stresses, because the use of deviatoric stresses can obscure the physical interpretation of calculated results. Deviatoric stresses are full stresses with the mean pressure removed, as called for by the constitutive relation. Deviatoric stresses thus include all three normal stresses. The force-balance equations are consequently more complex than when resistive stresses are used, except in the special case of plane flow. The partitioning of full stresses into lithostatic and resistive parts also allows a clearer distinction between action (driving stress) and reactions (gradients in resistive stresses). The action derives from gravity and the reactions from forces applied at the margins of the studied section of a glacier.

Our development follows closely the approach of prior work apart from the later introduction of deviatoric stresses and that boundary conditions at the upper surface are used instead of unknown basal conditions. Furthermore, there is no need to neglect the effects of the T-term, or “bridging effects”, as they are termed here. All stresses are included in the numerical calculation scheme and, as a result, the method is applicable to all glaciers: valley and outlet glaciers, ice shelves and ice streams, as well as inland ice.

Force Balance in Horizontal Directions

Consider the full stresses σij (i,j = x, y, z) acting on an element in the glacier with coordinate axes horizontal (x,y) and vertical (z). Balance of forces is expressed by:


in which g represents acceleration due to gravity, and ρ ice density. Inertial forces are neglected.

Full stresses are partitioned into lithostatic (L) and resistive components (R ij):


where δ¡j is the Kronecker delta and ft represents the elevation of the upper surface of the glacier. R¡j represents elements of the resistive stress tensor and their net effect is to oppose the motion of the glacier.

MacAyeal (1987) used a very similar partition. He worked with quantities integrated through the thickness measured perpendicularly to the bed. His “dynamic drag” corresponds to our horizontal resistive stress, and his “form drag” to the lithostatic stress.

Separating full stresses into resistive and lithostatic components, the equations of horizontal force balance become:


Integrating from the base of a column (z = b) to the surface (z = h), one obtains:


The integration can be taken from any general depth to the ice surface, but as a start, force balance is calculated for the full thickness.

The order of integration and differentiation is now reversed and basal drag and driving stress are defined. Basal drag, τbi (i = x,y) includes all basal resistances:


and the driving stress is the familiar formula:


where xx = x and xy = y (Note that the driving stress does not act on a specific surface and so it is not a true stress; rather it is the net horizontal effect of gravity per unit area in map view.) One then obtains new balance equations from Equations (7) and (8):


The upper boundary conditions require a stress-free surface, which in terms of resistive stresses is (i = x,y):


This greatly simplifies the balance equations to:


Equations (14) and (15) are the basic equations needed for calculating basal drag from driving stress and resistive stresses acting on vertical surfaces. The first term represents the contrast in resistive forces acting on the x-faces of an ice column, the second term is the imbalance in resistive forces acting on the y-faces. The third term describes the action through the driving stress, and the last term is friction at the base of the glacier.

The balance equations as derived here have close affinities with the results of earlier developments (for example, Raymond, 1980, equation (25)). Because no approximations are made in the above derivation, the balance Equations (14) and (15) are exact.

Bridging Effects

Because ice is very soft, the vertical normal stress is very nearly lithostatic and the vertical resistive stress, Rzz , is nearly zero. For simple calculations setting it to zero is appropriate, but this approximation need not be made in numerical computations.

Differences from lithostatic could be important at distances that are short compared to the ice thickness, or with certain special features such as ice falls. At small scales, part of the glacier can act like a bridge in which the vertical normal stress at the underside of the bridge span in less than the weight of ice above, while that under the abutment exceeds the weight of material above. This is illustrated in Figure 1. The variations in vertical normal stress from lithostatic lead to shear-stress gradients, as noted in the figure.

Fig. 1. The bridging effect; vertical normal stresses are more compressive under the abutment and less compressive under the span. There are associated variations in shear stress

This topic has been discussed many times before by, among others, Budd (1970a), Hutter (1983, p. 265), Kamb and Echelmeyer (1986), Nye (1969), Whillans (1987), and Whillans and Johnsen (1983). These authors concur that bridging effects are important only on short-distance scales (a few ice thicknesses at most). Our results of force-budget calculations along the Byrd Station Strain Network (Van der Veen and Whillans, 1989; hereafter referred to as part II) confirm this conclusion.

Thus, for first-order calculations, bridging effects may be neglected and the vertical normal stress taken to be lithostatic, or equal to the weight of ice vertically above. In this approximation, the vertical resistive stress, Rzz = σzz - L is zero, and


Differentiating with respect to z gives:


Comparing Equation (17) with the third balance Equation(3) shows that neglecting bridging effects is equivalent to the assumption that horizontal gradients in the vertical shear stresses are small compared to pg. This is satisfied if

When the balance equations are solved numerically, this approximation is not necessary, and the value of Rzz can be obtained from Equation (3) for vertical force balance. Separating the full stresses into resistive and lithostatic components, that equation becomes:


As before, the equation is vertically integrated, this time from a general depth z, to the surface, where the boundary condition is:


After reversing the order of integration and differentiation, the result is:


This expression, or the simplifying assumption Rzz = 0, is needed when force balance in terms of strain-rates is considered.

Relation between Resistive and Deviatoric Stresses

Constitutive relations generally link strain-rates to deviatoric stresses rather than to full stresses or resistive stresses. Deviatoric stresses are defined as full stresses minus the average normal stress (or spherical stress) which is supposed not to enter into the constitutive relation, except in a very minor way (for example, Hooke, 1981). Consequently, the trace of the deviatoric stress tensor is zero as it is for the strain-rate tensor (for incompressibility), and the two can be readily linked. Therefore, to apply a constitutive relation of this type, it is necessary to write resistive stresses in terms of deviatoric stresses, σ'¡j

The two schemes for separating full stresses are:


in which the spherical stress is defined as S = (σ xx + σyy + σzz)/3.

Eliminating σ¡j from these expressions:


Taking i = j = z


because, by the definition, the normal deviatoric stresses must sum to zero. Substituting this into the expressions (22) provides resistive stresses in terms of deviatoric stresses:


(Recall that, for simple calculations, it is a good approximation to set Rzz = 0.)

Replacing resistive stresses with deviatoric stresses according to equation (24), the equation for force balance in the x-direction becomes:


These two balance equations are exact and neither the x-nor y-axis need follow a flow line.

By imposing more or less severe restrictions, the balance equations reduce to those of earlier developments. Setting Rzz = 0 is equivalent to neglecting bridging effects and the last term on the left-hand sides (corresponding to the T-term) vanishes. For the case of σ'xy = 0 = σ 'yy Equation (25) reduces to that derived in Nye (1969). A valley “shape factor” is not included in equation (25) and (26); rather, the effect of side drag is included directly using side shear R or σ 'xy .

Application of Equation (25) or (26) need not be restricted to the flow of grounded ice. For instance, it is equally valid for ice shelves where basal drag is zero. Bridging effects can probably be neglected, and for an ice shelf in a parallel-sided bay, σ'yy = 0, and Equation (25) reduces to:


where the overbar denotes the vertical mean value, and H = h - b. In this, side-drag gradients are written as (Thomas, 1973a):

according to an approximation in which side drag is taken to vary linearly and W represents the half-width of the shelf. Equation (27) was used by Van der Veen (1987) in a numerical ice-shelf model. For interpretation of ice-shelf data, Thomas (1973b) used an x-integrated version of Equation (27) with sea-water pressure as a boundary condition.


The constitutive relation provides the linkage between deviatoric stresses and strain-rates. The most commonly used relation is due to Glen and Nye and reads (Paterson, 1981, p. 30–31):


in which the exponent n is usually taken as equal to 3, and the stiffness parameter B is very sensitive to temperature. The effective strain-rate is:


This non-linearity in the constitutive relation is the principal reason why numerical techniques are best used.

The horizontal resistive stresses can be expressed in terms of strain-pates using the relations (24) and the constitutive relation (28). This gives:


In terms of strain-rates, force balance for the horizontal directions (Equations (25) and (26)) is thus expressed by:


These equations are used in the numerical solution scheme discussed below. If integrations in Equations (31) and (32) are carried out from a general depth z to the surface, rather than over the entire thickness, the balance equations allow calculation of strain-rates at all depths without making any ad hoc assumption about their vertical variation. The requirement for this calculation scheme is that surface velocities be measured, as well as the glacier’s geometry.

For simple, quick calculations simplifying assumptions can be made. Rzz can be set to zero and strain-rates can be held constant with depth at their surface values. This simple calculation is compared with a fuller calculation for Byrd Glacier in part III of this series (Whillans and others, 1989). Although deep velocities cannot be calculated this way, the method yields deep stresses that are very similar to those obtained from the more exact calculation.

An account of the finite-difference solution scheme is given next, but in short, calculations start at the surface and proceed gradually downward (Fig. 2). At each depth, velocities at that depth and at all shallower strata are used to solve the balance equations for the shear strain-rates έxz and έ yz . These are used to calculate velocities at the next lower layer, and the calculations continue downward until the bed is reached.

Fig. 2. Illustration of the solution scheme for calculating velocities and stresses at depth. Panel A shows the input data (surface slope, ice thickness, surface velocities) and geometry of the model. From the velocities, longitudinal strain-rates are calculated (panel B) and, invoking the surface-boundary conditions, shear stresses (panel C). From these, the vertical shear in velocity follows, and velocities at the next depth layer can be calculated (panel D). The vertical velocity follows from the stretching in panel B and incompressibility. The next steps are to calculate longitudinal strain-rates (panel E) and to solve the force-balance equation for shear stresses (panel F). This procedure is repeated for all depth layers (panel G). Finally, when the bed is reached, basal drag and stress normal to the bed are calculated (panel H)

Vertical Scaling for Numerical Calculations

The equations expressing balance of forces (Equations (31), (32), and (20)) are very non-linear and cannot be solved analytically without simplifications. A scheme for numerical solution is, however, readily constructed.

For clarity, plane flow in the (x,z)-plane is considered in the description, with the x-axis horizontal and along the mean direction of flow, and z vertically upwards. This simplification implies that strain-rate components involving the cross-flow, y-direction are zero, and also that force budget in this direction need not be considered. Extension of the methods developed here to include the second horizontal direction is straightforward, and in part III the results of a full three-dimensional calculation are presented.

The first step toward solving the balance equations numerically is to define a dimensionless coordinate to account for variations in thickness H = h - b, along the flow line;


At the surface, z = h and s = 0, whereas at the bed, z = h - H and j = 1. Using this vertical coordinate, the model domain becomes a rectangular two-dimensional grid.

In the (x,s) coordinate system, Equation (1) for force balance in the horizontal direction becomes



A derivation of the relation between horizontal gradients in the (x,z) coordinate system, and those in the (x,s) coordinate system, can be found in Holton (1979, p. 21–22). Substituting the partitioning of full stresses into resistive and lithostatic components, and multiplying by H gives


Integrating this expression with respect to s, from a general depth, s, in the ice to the surface (s = 0) yields an equation for the shear stress at depth:


where the surface-boundary condition (13) has been used. This corresponds to Equation (14) and describes the balance of forces on a column from the surface to the level s.

To relate resistive stresses to observable velocities, the constitutive relation is used. From Equation (30) we find




Inserting these expressions into the balance Equation (36), force balance written in terms of strain-rates becomes


This is the basic equation used here. It allows calculation of the shear strain-rate, έχz , and hence the velocity variation with depth, for given driving stress, τdx , and longitudinal strain-rate, έχχ .

To calculate Rzz , Equation (20) describing force balance in the vertical direction is needed. Transforming this equation to the (x,s) coordinate system, and using Equation (38) to express Rxz in terms of strain-rates, the vertical balance equation becomes


Strain-rates are defined for the (x,z) coordinate system and these must be linked to the velocities which are computed in the (x,s) system. (Note, however, that the calculated velocities u and w are still parallel to the x- and z-axes, respectively; the coordinate transformation has not affected the directions of velocities, only the levels at which they are calculated.) This linkage to the (x,s) system makes the expression for strain-rate in the model domain somewhat more complex, but the benefits of using the (x,s) coordinate system are considerable.

The horizontal stretching rate in the (x,z) coordinate system is linked to velocity gradients through its definition:

and considering that the differential, du, can be written as

the horizontal stretching rate becomes:

The value of [∂s/∂x]z is obtained from the definition of s and using it gives


This cannot be solved directly because [∂u/∂s]x is related to έχz, and that is unknown until the balance Equation (39) is solved.

The relation between [∂u/∂s\x and ixz is obtained from the definition of the shear strain-rate:


In terms of gradients in the (x,s) coordinate system, the first term in the large brackets is:

The second term can often be neglected, especially when surface slopes are small, but that is not necessary. Following a development similar to that for έxx above, one obtains:

In this expression, the vertical derivative [∂w/∂s]x is linked to vertical strain-rate, and hence, because ice is incompressible, to horizontal strain-rate:


From Equation (42) we find


Substituting Equation (44) into Equation (41), the longitudinal strain-rate becomes:


This relation is used to solve the balance Equation (39). It relates the stretching έxx to the velocity gradients [∂u/∂x] and [∂w/∂x]s which are taken along surfaces of s = constant, and to the shear strain-rate έxz . In a similar fashion, the vertical resistive stress Rzz can be expressed in terms of these velocity gradients and ζxz . Thus, provided velocities at some depth s and above are known, the balance Equation (39) contains just one unknown, the shearing έxz, and therefore it is solvable by means of iteration, as discussed in the next section.

Numerical Solution

The focal point of the numerical scheme is solving the balance Equation (39) for the shear strain-rate έxz at general depth, s. As noted above, the solution can only be found after velocities at that depth layer, and at layers above, have been calculated. This means that calculations must start at the surface, where measured velocities are available, and proceed gradually downward, until the bed is reached.

The geometry of the numerical model is shown in Figure 2, panel A. Horizontal grid points are denoted by /, and depth layers by j, with j = 1 at the surface (where s = 0) and increasing downward to a maximum of j = J. The equations derived in the preceding section are readily transformed to this numerical grid using standard finite-difference methods (for example, Smith (1985), which also contains a discussion of round-off errors involved in this transformation).

Strain-rates at the surface are obtained from velocity gradients measured along the j = 1 (or s = 0) surface. Using the constitutive relation, the surface boundary condition (13) can also be written as


and combining with Equation (45) gives the horizontal stretching:


After calculating έ xx (ί,1), vertical shearing follows from Equation (46) and the vertical strain-rate from incompressibility:

Velocity change with depth for the surface stratum can now be calculated from Equation (44)

and from Equation (43)

Approximating the derivatives by forward differences, velocities at the second layer (j = 2) are:


where Δz(i) = H(i)/(J - 1) represents the vertical spacing between two adjacent layers.

Although έxz does not vary very much with depth near the surface, using Equation (48) can result in an over-or under-estimation of velocities at the second depth layer. Therefore, a two-step loop is included in the present model. First, all strain-rates at level j = 2 are calculated using the balance Equation (39), together with the velocities following from Equation (48). Next, improved values for these velocities are computed from central differencing


and these values are used to repeat the stress calculations at level j = 2.

At deeper layers (j ≥ 3) velocities can be calculated using central differencing:


However, applying Equation (50) to calculate u(i,j) can lead to “leap-frogging” velocities, with large values at odd levels and small values at even levels, or vice versa. Therefore, velocities are calculated using central differencing involving two higher layers.

The vertical gradient at the intermediate level, j - 3/2, is given by

and also by

Setting the right-hand sides equal:


Eliminating u(i,j - 2) from Equations (50) and (51) gives:


A similar expression is used to calculate the vertical velocity w(i,j)

The next step is to solve the balance Equation (39). Writing this equation for brevity as


where F represents the right-hand side of Equation (39) divided by A simple iteration scheme is used to find the solution, and this can be represented by:


in which superscripts denote the iteration steps. As an initial guess, the value of the layer above is used:


The solution scheme is summarized in Figure 2. Because the velocities at level j are known (panel D), longitudinal stretching (panel E) can be calculated from Equation (45), using Equation (55) for the shearing term. Similarly, Rzz can be expressed in terms of known velocity gradients and the unknown έxz, using Equation (40). Substituting these estimates in the right-hand side of the balance Equation (53) yields an improved value for έx and Rzz . Although the iteration is somewhat complicated because the balance Equation (53) also contains horizontal derivatives, convergence occurs rapidly. Taking as the convergence criterion that the relative difference between two subsequent calculations of εxz must be less than 0.1% along the entire flow line, the number of iteration steps is typically two near the surface, and five near the bed (where έxz varies more strongly with depth).

One may note that the calculation scheme can be inverted to start at the base and progressively work upward to the surface. Basal velocities must be specified, as well as the geometry of the section of glacier. Calculated results are basal drag and surface velocities. This procedure is described by Van der Veen (in press).

Concluding Remarks

The equations for force balance can be applied to the evaluation of the force budget of any glacier if the constitutive relation is prescribed. The scheme includes effects of lateral and longitudinal convergence or divergence of ice flow and of varying side drag along the glacier. The scheme also applies regardless of the azimuth of the horizontal coordinate axes, and thus is very useful for realistic situations in which the ice flow changes direction. Furthermore, the equations apply to floating ice (zero basal drag) as well as to grounded ice, and so a single set of equations can be used for glaciers that are partially afloat.

The calculation scheme presented here is straightforward and free from restrictive mathematical assumptions. Provided surface slopes, horizontal velocities at the ice surface, and ice thicknesses are known, stresses at all depths are calculated explicitly, without making a priori assumptions about velocity or stress variation with depth. In its finite-difference formulation, the solution scheme is very fast. A calculation on a horizontal grid with 47 grid points and 101 depth layers (for the Byrd Station Strain Network, in part II) requires about 15 min on an early model IBM-AT personal computer.

The primary motivation for this program is the availability of repeat photogrammetry of glaciers and the possibility of measuring hundreds of surface velocities distributed over a glacier. Most notably, such measurements have been carried out for Columbia Glacier (Meier and others, 1985), Byrd Glacier (Brecher, 1986), Ice Stream B (Whillans and Bindschadler, 1988), and Jakobshavns Isbras (personal communication from H.H. Brecher). In part III (Whillans and others, 1989), force-budget calculations for full three-dimensional flow of Byrd Glacier are discussed.

The method of calculating deep stresses can also be applied to smaller strain grids. Part II (Van der Veen and Whillans, 1989) applies the numerical solution scheme developed here to the flow along the Byrd Station Strain Network and discusses how the results are affected by uncertainties in the input data.


This work was supported by U.S. National Science Foundation grant DPP–83117235 and a Post-Doctoral Fellowship from The Ohio State University. We thank J. Bolzan, K. Hutter, and a referee for comments. Figures are by R. Tope and typing is by K. Doddroe and C Kinzelman. This is Byrd Polar Research Center contribution number 598.


Bindschadler, R. Stephenson, S.N. MacAyeal, D.R. Shabtaie, S.. 1987 Ice dynamics at the mouth of Ice Stream B, Antarctica. J. Geophys. Res., 92(B9), 88858894.
Brecher, H.H. 1986 Surface velocity determination on large polar glaciers by aerial photogrammetry. Ann. Glaciol., 8, 2226.
Budd, W.F. 1969 The dynamics of ice masses. ANARE Sei. Rep., Ser. A(IV). Glaciol. (Publ. 108.)
Budd, W.F. 1970 a. Ice flow over bedrock perturbations. J. Glaciol., 9(55), 2948.
Budd, W.F. 1970b The longitudinal stress and strain–rate gradients in ice masses. J. Glaciol., 9(55), 1927.
Collins, I.F. 1968 On the use of the equilibrium equations and flow law in relating the surface and bed topography of glaciers and ice sheets. J. Glaciol., 7(50), 199204.
Frolich, R.M. Mantripp, D.R. Vaughan, D.G. Doake, C.S.M.. 1987 Force balance of Rutford Ice Stream, Antarctica. International Association of Hydrological Sciences Publication 170 (Symposium at Vancouver 1987 — The Physical Basis of Ice Sheet Modelling), 323331.
Holton, J.R. 1979 An introduction to dynamic meteorology. Second edition. New York, Academic Press.
Hooke, R.Leb. 1981 Flow law for polycrystalline ice in glaciers: comparison of theoretical predictions, laboratory data, and field measurements. Rev. Geophys. Space Phys., 19(4), 664672.
Hutter, K. 1983 Theoretical glaciology; material science of ice and the mechanics of glaciers and ice sheets. Dordrecht, etc., D. Reidel Publishing Company.
Hutter, K. Legerer, F. Spring, U.. 1981 First–order stresses and deformations in glaciers and ice sheets. J. Glaciol., 27(96), 227270.
Kamb, B. Echelmeyer, K.A.. 1986 Stress–gradient coupling in glacier flow: I. Longitudinal averaging of the influence of ice thickness and surface slope. J. Glaciol., 32(111), 267284.
MacAyeal, D.R. 1987 Ice–shelf backpressure: form drag versus dynamic drag. In Van der Veen C.J. and J. Oerlemans, eds. Dynamics of the West Antarctic ice sheet. Proceedings of a Workshop held in Utrecht, May 6—8, 1985. Dordrecht, etc., D. Reidel Publishing Company, 141160.
MacAyeal, D.R. Bindschadler, R.A. Shabtaie, S. Stephenson, S. Bentley, C.R.. 1987 Force, mass, and energy budgets of the Crary Ice Rise complex, Antarctica. J. Glaciol., 33(114), 218230.
Meier, M.F. Rasmussen, L.A. Krimmel, R.M. Olsen, R.M. Frank, D.. 1985 Photogrammetric determination of surface altitude, terminus position, and ice velocity of Columbia Glacier, Alaska. U.S. Geol. Surv. Prof. Pap. 1258–F.
Nye, J.F. 1957 The distribution of stress and velocity in glaciers and ice–sheets. Proc R. Soc London, Ser. A, 239(1216), 113133.
Nye, J.F. 1969 The effect of longitudinal stress on the shear stress at the base of an ice sheet. J. Glaciol., 8(53), 207213.
Paterson, W.S.B. 1981 The physics of glaciers. Second edition. Oxford, etc., Pergamon Press.
Raymond, C.K. 1980 Temperate valley glaciers. In Colbeck S.C., ed. Dynamics of snow and ice masses. New York, Academic Press, 79139.
Robin, G. de Q. 1967 Surface topography of ice sheets. Nature, 215(5105), 10291032.
Smith, G.D. 1985 Numerical solution of partial differential equations: finite difference methods. Third edition. Oxford, Clarendon Press.
Thomas, R.H. 1973a. The creep of ice shelves: interpretation of observed behaviour. J. Glaciol., 12(64), 5570.
Thomas, R.H. 1973b. The creep of ice shelves: theory. J. Glaciol., 12(64), 4553.
Veen, C.J. 1987 Longitudinal stresses and basal sliding: a comparative study. In Van der Veen C.J. and J. Oerlemans, eds. Dynamics of the West Antarctic ice sheet: Proceedings of a Workshop held in Utrecht, May 6–8, 1985. Dordrecht, etc., D. Reidel Publishing Company, 223248.
Van der Veen, C.J. In press. A numerical scheme for calculating stresses and strain rates in glaciers. Math. Geol.
Van der Veen, C.J. Whillans, I.M.. 1989 Force budget: II. Application to two–dimensional flow along the Byrd Station Strain Network, Antarctica. J. Glaciol., 35(119), 6167.
Whillans, I.M. Bindschadler, R.A.. 1988 Mass balance of Ice Stream B, West Antarctica. Ann. Glaciol., 187193.
Whillans, I.M. Jezek, K.C.. 1987 Folding in the Greenland ice sheet. J. Geophys. Res., 92(B1), 485493.
Whillans, I.M. Johnsen, S.J.. 1983 Longitudinal variations in glacial flow: theory and test using data from the Byrd Station Strain Network, Antarctica. J. Glaciol., 29(101), 7897.
Whillans, I.M. Chen, Y.H. van der Veen, C.J. Hughes, T.. 1989 Force budget: III. Application to three–dimensional flow in Byrd Glacier, Antarctica. J. Glaciol. 35(119), 6880.