Introduction

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 (Reference ZienciewiczZienciewicz, 1971; Reference HuebnerHuebner, 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₋₅ (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

(1)

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

(2)

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 (Reference JohnsonJohnson, 1960; Reference HermannHerrmann, 1965).

For the present application Γ is defined as

(3)

where

(4)

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

(5)

where

(6)

which is equivalent to Glen’s law in the form

(7)

where

(8)

is the effective shear stress.

The temperature dependence of B is incorporated using a relation of the form (Reference LliboutryLliboutry, 1968) where . In these relations, ε̇_–5 is the strain-rate at a temperature of – 5°C, B_–5 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 Reference Mellor and TestaMellor 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 J_v in Equations (1) ; setting the derivative of J_v 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 Reference Nagtegaal, Nagtegaal, Parks and RiceNagtegaal 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 B_–5 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 (Reference NyeNye, 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 (Reference HookeHooke, 1973[a], Reference Hooke[b], Reference Hooke1976[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/m³, compared with 920 kg/m³ 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 (Reference HookeHooke, 1973[a], Reference Hooke[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 (Reference HookeHooke, 1973[b]) (see Table 1, column (2), and Fig. 1). Temperature measurements in the bore holes established the temperature distribution in the domain (Reference HookeHooke, 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 (Reference HookeHooke, 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 a^1/n (0.32 MN m⁻² a^1/n ) for the blue ice and 1.8 bar a^1/n (0.18 MN m^–2 a^1/n ) for the white.Footnote * 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 ₋₅ obtained for ice beneath the white ice band, 3.0 bar a^1/n (0.30 MN m⁻² a^1/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 a^1/n (0.44 MN m⁻² a^1/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 (Reference JonesJones, 1972), and surface velocities and strain-rates are known from repeated surveys of a series of 13 overlapping strain nets (Reference HoldsworthHoldsworth, 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. Footnote *

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 Reference HoldsworthHoldsworth’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 Reference Budd, Budd, Jenssen and RadokBudd 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 (Reference HoldsworthHoldsworth, 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 ₋₅ (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 (Reference BuddBudd, 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 ₋₅ might thus be higher. After a few trials, the distribution of B ₋₅ 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 Reference HoldsworthHoldsworth (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 (Reference HoldsworthHoldsworth, 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. y_b and y_s are the heights of the glacier bed and surface, respectively, above an arbitrary datum, and u and u_s are the horizontal velocities at heights y and y_s.

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 (Reference HookeHooke, 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,Footnote * 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 (Reference HuebnerHuebner, 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 Reference HoodHood (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 (Reference HermannHeinrich 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 (Reference HermannHeinrich 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 Reference RobinRobin (1955) and Reference BuddBudd (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⁻¹ 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⁻¹ 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 Reference HoldsworthHoldsworth’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⁻¹ at the distal edge to 0.019 deg m⁻¹ 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 Reference Budd, Budd, Jenssen and RadokBudd 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 programFootnote ^* 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 (Reference HoldsworthHoldsworth, 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 (Reference HookeHooke, 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 Reference HoldsworthHoldsworth’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 (Reference NyeNye, 1960; Reference Rasmussen and CampbellRasmussen and Campbell, 1973; Reference BuddBudd and Jenssen, [1975]; Reference Budd and JenssenBudd, 1975; Reference JenssenJenssen, 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.