Skip to main content Accessibility help


  • Access
  • Cited by 2



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

        Calculations of Velocity and Temperature in a Polar Glacier using the Finite-Element Method
        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.

        Calculations of Velocity and Temperature in a Polar Glacier using the Finite-Element Method
        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.

        Calculations of Velocity and Temperature in a Polar Glacier using the Finite-Element Method
        Available formats
Export citation


Numerical methods based on quadrilateral finite elements have been developed for calculating distributions of velocity and temperature in polar ice sheets in which horizontal gradients transverse to the flow direction are negligible. The calculation of the velocity field is based on a variational principle equivalent to the differential equations governing incompressible creeping flow. Glen’s flow law relating effective strain-rate ε̇ and shear stress τ by ε̇ = (τ/B) n is assumed, with the flow law parameter B varying from element to element depending on temperature and structure. As boundary conditions, stress may be specified on part of the boundary, in practice usually the upper free surface, and velocity on the rest. For calculation of the steady-state temperature distribution we use Galerkin’s method to develop an integral condition from the differential equations. The calculation includes all contributions from vertical and horizontal conduction and advection and from internal heat generation. Imposed boundary conditions are the temperature distribution on the upper surface and the heat flux elsewhere

For certain simple geometries, the flow calculation has been tested against the analytical solution of Nye (1957), and the temperature calculation against analytical solutions of Robin (1955) and Budd (1969), with excellent results.

The programs have been used to calculate velocity and temperature distributions in parts of the Barnes Ice Cap where extensive surface and bore-hole surveys provide information on actual values. The predicted velocities are in good agreement with measured velocities if the flow-law parameter B is assumed to decrease down-glacier from the divide to a point about 2 km above the equilibrium line, and then remain constant nearly to the margin. These variations are consistent with observed and inferred changes in fabric from fine ice with random c-axis orientations to coarser ice with single- or multiple-maximum fabrics. In the wedge of fine-grained deformed superimposed ice at the margin, B increases again.

Calculated and measured temperature distributions do not agree well if measured velocities and surface temperatures are used in the model. The measured temperature profiles apparently reflect a recent climatic warming which is not incorporated into the finite-element model. These profiles also appear to be adjusted to a vertical velocity distribution which is more consistent with that required for a steady-state profile than the present vertical velocity distribution.


Calculations of flow-rates and temperatures in glaciers and ice sheets can be done analytically only in the simplest of circumstances where simplifying assumptions about the geometry of the ice, the pattern of strain, the distribution of internal heat sources from straining, and the mechanical and thermal properties of the ice can be made. In real ice masses of complicated geometry and boundary conditions, numerical methods are necessary. In this paper we present a method, based on finite elements (Zienciewicz, 1971; Huebner, 1975), which represents a step toward treatment of the complexity of real ice masses. The two main restrictions are that the flow be planar, and the ice isotropic and quasi-viscous. Even with these restrictions, this is the most rigorous and practical method available for testing hypotheses about the flow behavior of glaciers in which complex geometries and spatial variations in rheological properties result in longitudinal variations in velocity and stress. Because of a lack of firm knowledge of the creep behavior of ice, no numerical model can give a priori accurate predictions of flow-rates and associated temperature distributions. However if, by judicious choice of flow-law parameters, a model such as the present one can be “tuned” to match field data, the data can then be extended to give realistic predictions of unmeasured quantities, and the sensitivity to changes in boundary conditions can be assessed.

In the present paper, our primary purpose is to give a brief introduction to the finite-element methods we use to calculate velocity and steady-state temperature distributions. As an example, we illustrate the method by applying it to parts of the Barnes Ice Cap where there is substantial geometrical and structural complexity (Figs 1 and 2).

Fig. 1. Cross-section of margin of Barnes Ice Cap along pole line B showing quadrilateral net used for finite-element calculations within this domain. Shading represents white ice. Large dots are nodes where velocity is known. Heavy vertical lines are bore holes. Variation in B−5 (value of B at – 5°C) used in finite-element calculation is shown at bottom.

Fig. 2. Cross-section along the trilateration net flow line, Barnes Ice Cap. See Figure 1 for explanation.

Calculation of Flow

The numerical method for the boundary conditions considered here is based on a functional


where A is the area of the longitudinal section in which planar deformation is to be determined, ρ is density, g is gravitational acceleration, u is the velocity vector, p is the mean compressive stress (pressure), and Γ is a function of the deformation-rate components

such that the components of deviatoric stress τij are given by


Equations (1) is appropriate to problems with stress-free conditions on part of the boundary (here the upper surface) and specified velocity on the remainder of the boundary (here no slip on the bed and specified velocity on any vertical edge); the solution is given by that velocity and pressure distribution which minimizes J v, and this is found by requiring δJ V = 0. This variational principle has been shown to be equivalent to the usual differential equations governing slow, incompressible, creeping flow (Johnson, 1960; Herrmann, 1965).

For the present application Γ is defined as




is the effective strain-rate and B and n are flow-law parameters. Summation over repeated subscripts is understood. Substitution of Equations (3) into Equations (2) gives




which is equivalent to Glen’s law in the form




is the effective shear stress.

The temperature dependence of B is incorporated using a relation of the form (Lliboutry, 1968) where . In these relations, ε̇5 is the strain-rate at a temperature of – 5°C, B5 is the flow-law parameter appropriate for that temperature, θ is the temperature in degrees Celsius (a negative number), and k is a constant with the values 0.34 for θ ⩾ –5°C and 0.15 for θ ⩽ – 5°C. These values of k are chosen to approximate the non-linear increase in ln ε̇ with decrease in reciprocal absolute temperature (1/T) observed by Mellor and Testa (1969, p. 135). A still lower value of k should be used if much of the glacier is colder than –15°C.

Quadrilateral finite elements are used to express the variational equations in a form which is practical for numerical solution. Within each quadrilateral element the distributions of velocity and pressure are defined in terms of the pressure parameters and in terms of the velocity components at the four corner nodes, using specific interpolating rules. The set of nodal velocity components and pressure parameters represent the discretized unknowns of the problem. These are chosen to minimize Jv in Equations (1) ; setting the derivative of Jv with respect to each unknown to zero gives a set of equations which can be cast in matrix form.

In the method used here, p is assumed constant in each quadrilateral. For evaluation of ε̇ (a Γand η and which are functions of ε̇) quadrilaterals are partitioned into two standard constant-strain-rate triangles along one diagonal allowing easy evaluation of the integral over the area of the element. Results for the two alternative partitionings are averaged. For evaluation of p ∇ · u and its integral over the area of an element, ∇ · u was taken to be equal to the average on the element. Although the evaluation of ε̇ and ∇ · u by different interpolation schemes lacks internal consistency and complete rigor, it has been found to give results essentially identical to a more rigorous treatment of a quadrilateral clement implemented in later programs in which the element is partitioned into four distinct constant strain-rate sub-triangles along the two diagonals, and the velocity at the diagonal intersection is chosen to give uniform dilatation rate, using a scheme developed by Nagtegaal and others (1974).

The program has three Fortran-coded modules: input, iteration, and graphical output, which communicate data through a random-access file. The input section does some internal reorganizing of node and element structure to minimize band-width of the matrix equations. The iteration section uses the strain-rate distribution associated with a trial solution to determine the effective viscosity distribution and then uses this viscosity distribution to generate an improved solution either by direct successive approximation or by Newton’s method. A modified Cholesky method is used to solve the matrix equations in any iteration step. The graphical section creates plots which show input data and calculated velocities and pressures throughout the longitudinal section of the glacier under consideration.

The input required for the program includes coordinates of nodes, M sequences of four-node numbers defining M quadrilateral elements filling the solution region, the temperature at each node, the specified velocity or stress at appropriate boundary nodes, and the flow-law parameters B5 and n. An initial velocity solution at interior nodes can also be given, and this can reduce the number of iterations required for satisfactory convergence. The program calculates u and v velocities at each node (u horizontal) and pressures in each element.

On the University of Minnesota Cyber 74 the Barnes Ice Cap Trilateration Net model (Fig. 2) with 250 nodes, 201 elements, and 51 fixed-velocity boundary nodes, takes 5.5 s per iteration. Convergence is usually satisfactory after six iterations, even if starting with zero initial values for the velocities. The program uses a total of approximately 121000 octal (= 41472 decimal) storage words (program, data, and temporary storage) for the model.

Tests of Flow Program

The flow program was tested in several ways. The behavior under various sets of simple boundary conditions and body force which theoretically give simple shear strain-rate, pure shear strain-rate, and hydrostatic stress, were all accurately reproduced by the calculated results. “Laminar flow” in an isothermal parallel-sided slab on a slope was also tested. Specific boundary conditions were a stress-free upper surface, no slip on the base, and periodic boundary conditions requiring the velocity to be the same at equivalent depths at the up-and down-stream faces of the slab. For n = 3 the solution accuracy depended on the number of nodes used over the depth, with five nodes being adequate for accuracy to two significant figures.

The most complex known analytical solution for planar flow of a power-law material is that for an isothermal, uniformly extending (or compressing), parallel-sided slab on a slope (Nye, 1957). The program was checked against this analytical solution for several combinations of slope, longitudinal strain-rate, and length-to-depth ratio of the solution region. In these tests, the boundary conditions were a stress-free upper surface and a specified velocity distribution on the ends and bottom of the slab. The boundary velocity values were determined from the analytic solution using the assumed longitudinal strain-rate. For n = 3, the calculated velocities were accurate to two or three significant figures when six or seven nodes were spaced over the depth of the slab.

Application of Flow Program

Two additional tests of the program were made using data from studies on the Barnes Ice Cap.

In the first of these we modeled the flow at the distal end of pole line B (Hooke, 1973[a], [b], 1976[b]) from bore hole B2 (Fig. 1) to the margin, a distance of 239 m. An interesting complication in the structure of this wedge-shaped domain is the presence of a band of white ice which owes its color to a high concentration of air bubbles. These bubbles reduce the density of the ice to 870 kg/m3, compared with 920 kg/m3 for the surrounding blue ice, and this results in a substantial reduction in the effective viscosity. Fabric studies have shown that this white ice has a strong single-maximum fabric. In contrast, the surrounding blue ice may have either a single- or a multiple-maximum fabric, except near the distal end of the wedge where the ice is much finer-grained and has a random fabric (Hooke, 1973[a], [b]).

This domain was divided into 96 elements with 116 nodes in such a way that the white ice (Fig. 1, stippled band) was distinct from the blue. Eighteen of the nodes within the domain and five nodes on its up-glacier boundary were placed to coincide with points where the velocity is known from triangulation and bore-hole deformation data (Hooke, 1973[b]) (see Table 1, column (2), and Fig. 1). Temperature measurements in the bore holes established the temperature distribution in the domain (Hooke, 1973[a], fig. 5). A stress-free boundary condition was used on the upper surface and a no-slip boundary condition on the frozen bed.

Table 1. Comparison of measured and calculated velocities

Analysis of the bore-hole deformation measurements and of deformation data from a series of 67 strain nets in the wall of a 125 m tunnel yielded estimates of the values of the flow law parameters n and B (Hooke, 1973[b]). n was estimated to be 1.65 ± 0.12. The most reliable values of B were thought to be those resulting from the analysis of bore holes B1 and B2. The resulting values, adjusted to a temperature of –5°C, are 3.2 bar a1/n (0.32 MN m−2 a1/n ) for the blue ice and 1.8 bar a1/n (0.18 MN m–2 a1/n ) for the white.* Use of these values in the finite-element program did not give good agreement between measured and calculated velocities, however (Table 1, column (3)) ; calculated velocities near the margin were too high. This suggested that some higher values of B obtained from analysis of bore holes B0 and B1x might, in fact, be real. It was therefore assumed that the value of B 5 obtained for ice beneath the white ice band, 3.0 bar a1/n (0.30 MN m−2 a1/n ), should apply midway between bore holes B1 and B2, and that the value obtained for ice near the bottoms of bore holes B0 and B1x, 4.4 bar a1/n (0.44 MN m−2 a1/n /), should apply midway between these two holes. The increase in B 5 was assumed to be linear as shown in Figure 1. This down-glacier increase in B is plausible in view of the distribution of fabric discussed above. As can be seen from Table 1 (column (4)), use of these values of B in the finite element program resulted in somewhat better agreement between calculated and observed velocities. In many cases the difference is within the limits of uncertainty in the measurements. The largest discrepancies are near the down-glacier limit of the white ice band (nodes P3 and B1: 0.0) where complications resulting from the viscosity contrast may be significant.

Numbers in columns (3), (4), (7), and (8) are delta values obtained by subtracting the calculated velocity from the known velocity. Negative u velocities are in the down-glacier direction and negative v velocities are downward.

For our second test we modeled the Trilateration Net flow line from the divide to the margin. Bed topography along this 10 km flow line has been determined by radio depth-sounding (Jones, 1972), and surface velocities and strain-rates are known from repeated surveys of a series of 13 overlapping strain nets (Holdsworth, 1975). In addition, temperature measurements (unpublished data of R. LeB. Hooke) have been made in a series of six 25 to 30 m bore holes approximately evenly spaced along the flow line and in four deeper holes, T061, T081, T091, and T0975, which penetrated well into the white ice at the base of the glacier (Fig. 2). (Hole numbers and numbers of labeled nodes in Figure 2 give the approximate distance from the divide–061 being 6.1 km; 0975, 9.75 km; and so forth). Finally, study of cores from holes T061, T081, and T0975 indicates that the basal white ice with a single maximum fabric is overlain by a 40 to 60 m thick layer of ice with a multiple (2 to 4) maximum fabric, which in turn is overlain by a 75 m thick layer of ice with a broad single maximum fabric. Still higher in the glacier, crystals are smaller and c-axes are randomly oriented. *

This domain was divided into 201 elements with 251 nodes (Fig. 2). Velocities are known at 14 of the surface nodes (Table 1, column (6)) from Holdsworth’s (1975) surveys. Nodal temperatures were estimated from the measured temperature distribution in the down-glacier half of the domain, and from a temperature distribution calculated using Budd and others’ (1971) column model in the up-glacier half. A velocity boundary condition was used at the up-glacier end of the flow line with u = 0 everywhere, and with v = 0.29 m a–1 at the surface (Holdsworth, 1975) and decreasing linearly with depth. A no-slip boundary condition was used at the frozen base and a zero-stress boundary condition along the surface.

Because ice fabrics in bore holes T061 and T081 are comparable to those between holes B1 and B2 on pole line B, the initial finite-element calculation for this problem was done with the values of n (1.65) and B 5 (3.2 in the blue ice and 1.8 in the white) obtained from analysis of deformation of these latter holes, as described above. This resulted in horizontal velocities which were substantially too high. A higher value of n would be consistent with the higher stresses in this part of the glacier (Budd, 1969, p. 21), so a value of 2.6 was tried. This reduced the velocities, but the differences between observed and calculated velocities were not uniformly distributed (Table 1, column (7)). Because the largest discrepancies were in the up-glacier half of the domain, it was concluded that the ice here might not have recrystallized as extensively as that farther down glacier, and that B 5 might thus be higher. After a few trials, the distribution of B 5 shown in Figure 2 was found to give good agreement between calculated and observed horizontal velocities (Table 1, column (8), and Fig. 3).

Fig. 3. Horizontal velocity calculated from finite-element program using n = 2.6 and variation of B shown in Figure 2, compared with measured horizontal velocity. Small dots are results calculated at nodes on surface. Open circles show values measured on modeled flow line. Xs show average of values measured on flow lines approximately 0.5 km east and west of modeled flow line. Measured values from Holdsworth (1975).

Calculated and observed vertical velocities still did not agree well, however (Fig. 4), and several factors may contribute to these discrepancies. Relatively small longitudinal variations in horizontal velocity, and hence in longitudinal strain-rate, possibly arising from small calculation errors associated with the descretized finite element network, or possibly caused by imperfect adjustment of the assumed surface profile to the assumed bed profile, may lead to substantial variations in vertical velocity. This leads to an artificial high-frequency noise, which could be eliminated by rather small adjustments of the upper surface topography. Small transverse compressive strain-rates also exist (Holdsworth, 1975) and are probably responsible for the observed velocities being systematically higher than the calculated results (Fig. 4).

Fig. 4. Comparison of calculated and observed vertical velocity. Symbols are the same as in Figure 3.

A further comparison of calculated and measured velocity is possible in hole T061 (Fig. 5). The vertical variation of horizontal velocity shows a tendency to be less concave near the ice divide than distant from it. At distances more than 3 km from the divide, the profiles are fairly consistent, but are still less concave than observed at 6.1 km. This is suggestive of either a higher value of n or of a decrease in B with depth in the blue ice. This latter possibility is suggested by the depth variation of fabric described above. No attempts were made to adjust n or B to produce a better fit to this vertical profile.

Fig. 5. Dimensionless velocity profiles along trilateration net flow line. Points labeled 010 are from profile 1 km from divide and so forth (see text). Measured profile in hole T061 is given for comparison. yb and ys are the heights of the glacier bed and surface, respectively, above an arbitrary datum, and u and us are the horizontal velocities at heights y and ys.

Discussion of Flow Calculations

The main uncertainty in the flow calculations lies in the discretization error associated with the relatively small number of elements and their specific arrangement in the models. A larger number of elements, or possibly a more nearly optimal arrangement for the number used, is possible; the main limitations are only the effort in setting up input, which is done manually because of the complicated structure, and computer time, the expense of which climbs dramatically with increasing number of nodes. Based on the tests against simple analytic solutions, we believe the results of the fairly coarse grid we used should be accurate to better than 10% except possibly very near the margin where the number of grid points over depth drops. This is adequate to draw some conclusions.

We can conclude with some certainty that the flow-law parameters deduced for the blue and white ices by analysis of bore-hole and tunnel deformation data (Hooke, 1973[b]) do not apply homogeneously and uniformly to these ice types over the whole ice-sheet section. Based on a number of trials with various assumed constant values of B and n, we conclude with somewhat less certainty that no homogeneous and uniform power-type flow law can adequately explain the observed motion. Therefore, we are led further to conclude that, in addition to the already documented differences between blue and white ice, the known structural heterogeneity and non-uniformity of the blue ice is also reflected in flow behavior which varies from place to place. The values of n and B assumed in the models seem reasonable in view of what is known about the structure of the ice,* but we do not suggest that these values represent a unique or best interpretation. A depth variation also seems likely, and further effort is needed to include this in the calculations. In this regard, it is important to recognize a certain inconsistency in the efforts described herein; spatial variation of B is attributed to variation of c-axis fabric, which implies a plastic anisotropy not taken into account by the calculation procedure.

In addition to these specific results regarding the behavior of sections of the Barnes Ice Cap, the calculations suggest some possibly more general characteristics of the behavior of ice caps spreading from a ridge. The tendency for depth profiles of horizontal velocity to be less concave near the divide than on the flanks is probably related to the symmetry condition at the divide rather than to any specific details of the geometry or flow-law assumptions. The irregular variation of vertical velocity along the surface suggests a strong sensitivity of vertical velocity to small perturbations on the upper surface. If this is true, prediction of vertical strain-rate and the relation between age and depth by numerical methods must be done cautiously. These features of flow near a divide have been found in calculations using the geometry and temperature distribution of the Devon Island ice-cap divide, and also in calculations using an idealized divide geometry, and are probably quite general (unpublished calculations and manuscript in preparation by C. F. Raymond).

A disturbing result of the finite-element calculations is that pressures in the elements oscillate longitudinally. In the trilateration net section longitudinal oscillations of about 20% were calculated locally in a series of elements along the bed. An overall down-glacier decrease in pressure is consistent with the decreasing ice thickness, but the oscillation is an artifact of the calculation. Similar but smaller oscillations occurred in a few places in the solution for pole line B. In the solution for the slab on a slope they were almost non-existent, being present only in the basal layer of elements, and then with a relative amplitude of less than 1%.

Method for Temperature Calculations

The finite-element equations for description of the temperature distribution within the glacier are developed from the steady-state field equation with combined advection and diffusion:

where θ is the temperature, k the thermal conductivity, c v the heat capacity, and u the velocity. This assumes no heat production from chemical or phase-change effects.

Galerkin’s technique (Huebner, 1975, chapter 4) is used to develop the finite-element expressions for this governing differential equation. The Galerkin development introduces the natural boundary conditions to allow specification of temperatures or temperature gradients on any boundary. Inclusion of the advective terms causes matrices to be non-symmetric, but they remain banded. The frontal solution technique as developed for finite-element applications by Hood (1976) is used to solve the equations.

The heat generation dτ in each element was calculated from the effective strain-rate. The components of the strain-rate tensor were calculated from their definition in terms of the velocity derivatives, which in turn are given by the specified nodal velocities. From Equations (4), (5), (6), and (8) dτ = 2B ε̇(1+1/n). Because B is temperature dependent, an iterative solution method is required, but convergence is rapid.

Most of the input for the temperature program is virtually identical to that for the flow program. The main differences are that: (1) nodal velocities must be specified; (2) the boundary conditions are now either a temperature or a temperature gradient; and (3) two additional material properties, the thermal conductivity and the heat capacity of the ice, must be specified.

It has been recognized that numerical solution of the heat-flow equation can present serious difficulties when the advective term is important (Heinrich and others, 1976). These difficulties stem from the combination of the essentially elliptic and parabolic nature of the two terms, and result in oscillation in the solution whenever the mesh size exceeds a certain critical value. Such oscillations with amplitudes up to 0.5 deg in the longitudinal direction, occurred initially in the trilateration net model. To reduce this oscillation to a more acceptable level, up-winding was incorporated into the program (Heinrich and others, 1976).

On the University of Minnesota Cyber 74 the trilateration net model required 7.4 s per iteration. Convergence is usually satisfactory after three iterations. The program uses approximately 50000 octal (= 20480 decimal) storage words for this application.

Tests of Temperature Program

The temperature program was tested against analytical solutions of Robin (1955) and Budd (1969, p. 67, equation (26)). These tests were done on a slab 1000 m thick with a basal temperature gradient of 0.02 deg m–1. In one test, the horizontal velocity u was set to zero and v varied linearly from 0 at the bottom to 0.3 m a−1 downward at the surface, giving a temperature difference across the slab of 8.9 deg. The calculated profile was within 0.05 deg of Robin’s analytical solution. Slightly warmer temperatures were calculated by the finite-element program owing to mechanical generation of heat which is not taken into account by Robin’s solution. In a second test, u was set equal to 5 m a−1 on the surface and decreased to zero at the bed in accord with Glen’s law, and v was set to zero everywhere, giving a temperature difference across the slab of 25.0 deg. The resulting profile was within 0.04 deg of Budd’s analytical solution.

Application to Barnes Ice Cap

The temperature distribution along the trilateration net flow line (Fig. 2) was initially calculated using input velocities deduced from Holdsworth’s (1975) surface measurements. The variation of u with depth was determined from Glen’s law with n = 3 and with B adjusted to give u = 0 on the bed; v was assumed to decrease linearly with depth. The temperature on the upper boundary was measured in shallow bore holes, and the temperature gradient at the bed was allowed to vary from 0.015 deg m−1 at the distal edge to 0.019 deg m−1 at the divide in accord with measurements in holes T061 and T091.

The finite-element model gave a temperature distribution in the up-glacier half of this domain which was in reasonable agreement with Budd and others’ (1971) column model. Differences between the two models were generally less than 0.5 deg, and can be accounted for by slight differences in the vertical velocities assumed for the two models and by the approximations used in the column model to deal with horizontal advection. In the down-glacier part of the domain, however, good agreement between measured and calculated temperatures was not obtained. In T061 the calculated temperature profile was systematically about 1.4 deg warmer than the measured profile over the bottom 75% of the glacier, and in T091 the discrepancy was even greater (Fig. 6, long-dashed curves). The magnitude of the discrepancies suggested that the measured profiles were adjusted to either lower (i.e. more negative or downward, or less positive) vertical velocities, than used in the finite-element calculation, or to lower horizontal advection (lower u or dθ/dx).

Fig. 6. Measured and calculated temperature profiles in bore holes T061 and T091. Long-dashed curve was calculated using measured vertical velocities, whereas short-dashed curve was calculated using vertical velocities required for a balanced mass balance with present accumulation pattern, horizontal velocities, and surface profile (see text).

Modeling with the use of a one-dimensional finite-difference program * using observed values of u and dθ/dx, and consideration of the sensitivity of the vertical velocity to small changes in longitudinal strain-rate suggested that uncertainties in the vertical velocity were most likely. If the accumulation-rate (negative for ablation) ȧ is assumed to be given by ȧ = f (x) + c, where x is the longitudinal coordinate and f (x) is obtained from recent measurements (Holdsworth, 1975; and unpublished measurements of R. L. Hooke), the constant c can be adjusted to give a balanced budget. Then v = ȧ + u tan α, where α is the surface slope, yields a new vertical velocity distribution. Use of this velocity distribution considerably improved the agreement between calculation and measurement (Fig. 6, short-dashed curves). Discrepancies still exist, however, and further modeling will be necessary to identify possible causes of these discrepancies.

One factor not taken into consideration in the finite-element model is surface warming in response to climatic warming, possibly accentuated by the effects of percolating melt water (Hooke, 1976[b]). The reverse curvature near the top of the T061 profile (Fig. 6) implies such warming, and finite-difference modeling suggests that this reversal could have been caused by a warming of slightly less than a degree about a decade ago.

A tentative conclusion based on work done to date is that vertical velocities at the surface at present are higher than those to which the temperature distribution and surface profile of the ice cap are adjusted. This contrasts with Holdsworth’s (1975, p. 17) conclusion that “dynamically the ice cap seems to be reacting to a lower accumulation-rate and higher ablation rate than are currently recorded”. The temperature modeling suggests instead that dynamically the ice cap is in fact responding to a higher-than-average accumulation-rate sometime in the recent past.

Perspective on Methods

It is worthwhile contrasting the methods described here with other models of glaciers and ice sheets (Nye, 1960; Rasmussen and Campbell, 1973; Budd and Jenssen, [1975]; Budd, 1975; Jenssen, 1977). These latter models are time-dependent; they calculate the redistribution of mass with time as their principal goal, and as such are concerned more with ice discharge and less with details of the velocity distribution. With this goal it is appropriate to introduce averages over depth in cross-sections of a glacier and assumptions about vertical variation of velocity in order to arrive at a practical computational scheme which can be stepped through long intervals of time.

Except for the restrictions to planar deformation and isotropic quasi-viscous flow, the model we present avoids such simplifications and is limited in accuracy only by the density of grid points and knowledge of the correct flow law. It gives the velocity distribution in detail. Conceivably, it could be used for time-dependent models and the actual computer programming is not difficult, but this is not very practical because of the large amount of computer time required for each time step. The utility of the method lies more as a means to test specific hypotheses against measurements. For example, the methods we describe could be used to test the simplified time-dependent models mentioned above at any given time step.

Concluding Statement

An obvious conclusion that can be drawn from the modeling we have done so far is that in the absence of measurements which can be used to “tune” the models, the calculated velocities and temperatures must be used with considerable caution. Also, tuning such models by adjusting temperature gradients or spatial variation of flow-law parameters, for example, to better match field data, is a somewhat uncertain process. We have successfully improved the agreement between calculation and measurement in the flow models for the Barnes Ice Cap through several trials of B and n, and have arrived at plausible values of these constants with spatial variations which are reasonably compatible with observation and intuition. Better agreement could undoubtedly be obtained with further trials, but until experimental data on the variations of B and n with fabric and grain size are available, the range of possible values for these parameters is too large to make such an exercise productive. Better agreement between calculated and measured temperatures probably could also be obtained by varying the boundary conditions and the vertical velocities. In either case, however, the uniqueness of the solutions is questionable in the absence of field and laboratory support for any assumptions that are made.

One of the most important results of the modeling we have done so far is that attention has been sharply focused on the need for a more comprehensive constitutive relation between stress and strain-rate. This relation must include not only the effects of fabric and grain size on the · flow-rate under a given stress, but must also address the problem of the progressive development of fabrics with increasing total strain.


This research was supported by National Science Foundation grant DPP74-19075 to the University of Washington, grants GA-42728 and EAR77-21098 to the University of Minnesota, and by the University of Minnesota Computer Center. Logistical support for some of the field work was provided by the Glaciology Division, Environment, Canada.

page note * A bore hole, B22, drilled in 1974, after publication of the results of the original work in this area, established that the ice was only 31.3 m thick at the upper end of this section, rather than 35.5 m as originally assumed. The bore-hole deformation measurements have been re-analyzed for the present study, taking this new datum into consideration.

page note * Paper entitled “Changes in ice fabric and texture with depth and with distance along a flow line”, by R, L. Hooke and P. J. Hudleston, is in preparation.

page note * Since completion of this paper, fabric data were obtained from a 300 m hole at 020 (Fig. 2)A random fabric on the upper part of this hole gives way directly to a multiple-maximum (2-3 maxima) fabric at a depth of about 130 m. This is not consistent with the variation in B_5 assumed in the calculations.

page note * Paper entitled “Temperature profiles in the Barnes Ice Cap, Baffin Island, Canada, and heat flux from the subglacial terrain”, by R. L, Hooke, R. J. Gustafson, and E. C. Alexander, Jr.


Budd, W. F. 1969. The dynamics of ice masses. ANARE Scientific Reports. Ser. A(IV). Glaciology. Publication No. 108.
Budd, W. F. 1975. A first simple model for periodically self-surging glaciers.. Journal of Glaciology, Vol. 14, No. 70, p. 321.
Budd, W. F., and Jenssen, D. [1975.] Numerical modelling of glacier systems. [Union Géodésique et Géophysique Internationale. Association Internationale des Sciences Hydrologiques. Commission des Neiges et Glaces.] Symposium. Neiges et glaces. Actes du colloque de Moscow, août 1971, p. 257–91. (IAHS-AISH Publication No. 104.)
Budd, W. F., and others. 1971. Derived physical characteristics of the Antarctic ice sheet, by Budd, W. F., Jenssen, D., and Radok, U.. ANARE Interim Reports. Ser. A(IV). Glaciology. Publication No. 120.
Heinrich, J. C, and others. 1976. An “upwind” finite element scheme for two-dimensional convective transport equation, by Heinrich, J. C., Huyakorn, P. S., Zienkiewicz, O. C., and Mitchell, A. R.. International Journal for Numerical Methods in Engineering, Vol. 11, No. 1, p. 131–43,
Hermann, L. R. 1965. Elasticity equations for incompressible and nearly incompressible materials by a variational theorem. American Institute of Aeronautics and Astronautics Journal, Vol. 3, No. 10, p. 18961900.
Holdsworth, G. 1975. Deformation and flow of Barnes Ice Cap, Baffin Island. Ottawa, Environment Canada. Water Resources Branch. Inland Waters Directorate. (Scientific Series, No. 52.)
Hood, P. 1976. Frontal solution program for unsymmetric matrices. International Journal for Numerical Methods in Engineering, Vol. 10, No. 2, p. 370–99.
Hooke, R. L. 1973[a]. Flow near the margin of the Barnes Ice Cap and the development of ice-cored moraines. Geological Society of America. Bulletin, Vol. 84, No. 12, p. 3929–48.
Hooke, R. L. 1973[b]. Structure and flow in the margin of the Barnes Ice Cap, Baffin Island, N.W.T., Canada. Journal of Glaciology, Vol. 12, No. 66, p. 423–38.
Hooke, R. L. 1976[a]. Near-surface temperatures in the superimposed ice zone and lower part of the soaked zone of polar ice sheets. Journal of Glaciology, Vol. 16, No. 74, p. 302–04.
Hooke, R. L. 1976[b]. Pleistocene ice at the base of the Barnes Ice Cap, Baffin Island, N.W.T., Canada. Journal of Glaciology, Vol. 17, No. 75, p. 4959.
Huebner, K. H. 1975. The finite element method for engineers. New York, John Wiley and Sons, Inc.
Jenssen, D. 1977. A three-dimensional polar ice-sheet model. Journal of Glaciology, Vol. 18, No. 80, p. 373–89.
Johnson, M. W. 1960. Some variational theorems for non-Newtonian flow. Physics of Fluids, Vol. 3, No. 6, p. 871–78.
Jones, S. J. 1972. Radio depth-sounding on Meighen and Barnes Ice Caps, Arctic Canada. Ottawa, Environment Canada. Water Resources Branch. Inland Waters Directorate. (Scientific Scries, No. 25.)
Lliboutry, L. A. 1968. Steady-state temperatures at the bottom of ice sheets and computation of the bottom ice flow law from the surface profile. Journal of Glaciology, Vol. 7, No. 51, p. 363–76.
Mellor, M., and Testa, R. 1969. Effect of temperature on the creep of ice. Journal of Glaciology, Vol. 8, No. 52, p. 131–45.
Nagtegaal, J. C., and others. 1974. On numerically accurate finite element solutions in the fully plastic range, by Nagtegaal, J.C., Parks, D. M., and Rice, J. R.. Computer Methods in Applied Mechanics and Engineering, Vol. 4, No. 2, p. 153–77.
Nye, J. F. 1957. The distribution of stress and velocity in glaciers and ice-sheets. Proceedings of the Royal Society of London, Ser. A, Vol. 239, No. 1216, p. 113–33.
Nye, J. F. 1960. The response of glaciers and ice-sheets to seasonal and climatic changes. Proceedings of the Royal Society of London, Ser. A, Vol. 256, No. 1287, p. 559–84.
Rasmussen, L. A., and Campbell, W. J. 1973. Comparison of three contemporary flow laws in a three-dimensional time-dependent glacier model. Journal of Glaciology, Vol. 12, No. 66, p. 361–73.
Robin, G. de Q. 1955. Ice movement and temperature distribution in glaciers and ice sheets. Journal of Glaciology, Vol. 2, No. 18, p. 523–32.
Zienciewicz, O. C. 1971. The finite element method in engineering science . London, McGraw-Hill Book Co., Inc.


L. W. Morland: It is possible to prescribe accumulation/ablation-rates to conform with smooth velocities observed on the surface, and then determine the free surface, but this is a much more difficult computation.

C. F. Raymond: A smooth distribution of emergence velocity along the surface could be obtained by letting the surface adjust in time to come into balance with a prescribed, smoothly varying distribution of net balance. This could be done with our method, but it would involve stepping in time and would probably take too many time steps.

Morland: Zero velocity at the base is incompatible with finite ablation at the margin, so there is a singularity in the solution which is not described by the numerical method. Or, as in this case, the solution must imply zero ablation at the margin which is inconsistent with observation.

Raymond: We do not prescribe net balance. The free upper surface moves up and down according to the velocity determined by the geometry, assumed ice properties, and the idealization by finite elements. Right at the edge, the velocity is zero, and of course this geometry could not be maintained against ablation. Since we seek the velocity at a given time under a given geometry I do not think this is a serious problem. If we were to time iterate as envisaged in my answer to your first question, this would be a definite problem.

R. Brepson: I suppose that you compute the plastic deformation with the deviatoric stress at each step. Thus you are right to use high isoparametric elements. How many iterations are necessary to achieve convergence?

Raymond: Our actual iteration procedure involves computation of effective viscosity from effective strain-rate associated with a trial velocity distribution. This effective viscosity distribution is used to recalculate an improved velocity distribution. Once the trial solution is close enough to the real one (1 to 10%), then the Newton-Raphson method we use converges very rapidly. Only several iterations are required to achieve 10−6 to 10−8 fractional change from iteration to iteration.

Brepson: What is the coupling between the temperature solution and the deformation solution? Is it the dependence of viscosity on temperature?

Raymond: The coupling comes through the effect of temperature on effective viscosity, which affects the mass flow and velocity which causes advective transfer of heat.

O. Orheim: Your two measured temperature profiles seem remarkably similar to the depth measured of about 100 m, even though one was in the ablation area and one in the accumulation area, and we tend to think that the heat input into the latter area must be much greater because of refreezing of melt water.

N. W. Young : If you consider the full equation you are using, the inclusion of the second derivative in the horizontal ∂2θ/∂x 2 of the temperature admits oscillating solutions in the horizontal because of the relative proportions of the dimensions. Removal of that term, as is possible in finite-difference techniques, would give a more reasonable solution without that problem.

Raymond : The steady-state heat-flow equation with both diffusion and advection terms has a combined elliptic-hyperbolic character which is apparently the source of the problem. Intuitively I would expect that the term ∂2θ/∂x 2 would tend to even out oscillations by diffusion and therefore help the situation. If it is advantageous to leave it out (without loss of accuracy in the glacier calculations where horizontal conduction is negligible), this is more difficult to accomplish in the finite-element method than the finite-difference method.

D. A. Yuen: What is the maximum effective viscosity encountered in your calculations? Was any cut-off effective viscosity used to stabilize the Newton-Raphson scheme used in the set of non-linear algebraic equations ? What is the maximum difference in effective viscosity among all sets of nodal points?

Raymond: No effective viscosity cut-off was used to stabilize the Newton-Raphson scheme. Off-hand I do not know the maximum viscosity contrast in the Barnes Ice Cap solutions.