Introduction

Glaciers and other large ice masses provide detailed and continuous records of past climatic and atmospheric chemistry over a period of approximately 100000 years. For ice-core information to be useful, the age of the glacier ice must be established. There are four basic methods for determining the age of glacier ice: (1) radioisotope methods; (2) use of reference horizons; (3) counting annual layers; (4) use of flow models (Reference PatersonPaterson, 1981). Each of these methods has advantages and disadvantages that are associated with the respective techniques for dating. A review of these techniques is not intended; instead, the reader is referred to the aforementioned reference for details.

This paper outlines a procedure for determining the age of glacier ice based on predictions of the velocity field and associated particle paths by using velocity-pressure and streamline finite-element formulations. The finite-element models presented in this paper offer distinct advantages over existing flow models for determining particle paths, and thus for establishing the age of glacier ice. The advantages associated with the finite-element method are:

• (1) More realistic velocity fields can be determined to establish particle paths.

• (2) Variations in material properties and temperatures can be accommodated.

• (3) Glaciers with irregular geometry and complex boundary conditions can be modelled.

It is assumed that ice creep can be modelled by using contemporary continuum-mechanics concepts for planar flow and that the variation in ice-creep properties can be established through optimization procedures which estimate the flow parameters by minimizing the error between predicted and observed surface velocities. Of course, in order to predict particle paths based on steady-state flow fields, we must assume that the ice mass has been deforming under steady-state conditions over a period of time greater than the age of ice in which we are interested.

Examples are presented which clearly demonstrate the applicability of the proposed procedure to predict particle paths and isochrones.

Background Preliminaries

The general procedure for establishing the age of ice via flow models, involves:

• (1) Determining velocity field of ice mass.

• (2) Estimating the particle paths from the velocity field.

• (3) Integrating along particle paths to obtain the age of ice along the particle paths (see Fig. 1).

• (4) Joining points of equal time between particle paths to establish isochrones.

Fig. 1. Particle path for typical glacier.

Many of the past methods for establishing particle paths and isochrones have relied on simplifying assumptions regarding the flow field which often have ignored significant features of glacier flow. For example, Budd’s model (Reference Budd, Budd, Jenssen and RadokBudd and others, 1971) assumed glacier ice flows as a column with all shear concentrated at the bottom, resulting in almost uniform horizontal velocities and constant vertical strain-rate with depth. While such dynamics may reflect the flow behaviour of ice masses having significant basal sliding or high basal shearing due to warmer basal temperatures, such a flow field, in general, presents an over-simplification. For further information on dating glacier ice with simplified flow models. the reader is referred to Reference NyeNye (1963), Reference Dansgaard and JohnsenDansgaard and Johnsen (1969), Reference Philberth and FedererPhilberth and Federer (1971), and Reference Budd, Budd, Jenssen and RadokBudd and others (1971).

Recently, Reference Reeh, Reeh, Johnsen, Dahl-Jensen, Langway, Langway, Oeschger and DansgaardReeh and others (1985) modelled the Dye 3 area, Greenland, by combining finite-difference and perturbation techniques to establish flow-line and isochrone patterns. Although their approach is a departure from the simplified ones which have generally been adopted by glaciologists in the past, the authors believe that it is more complicated and not as versatile as the approach expounded in this paper.

Finite-Element Model

Owing to the advent of the finite-element method in the 1950s, many complex boundary-valued problems, for which closed-form analytical solutions do not exist, have been solved. The method has been successfully applied to many fields of endeavour, and more recently has been applied to the problem of glacier dynamics (Reference Hooke, Hooke, Raymond, Hotchkiss and GustafsonHooke and Others, 1979; Reference RaymondRaymond, 1983; Reference StolleStolle, unpublished). This section summarizes the finite-element models used in this study for the determination of particle paths.

To solve for particle paths, it is possible to approach the problem in terms of a stream-function formulation as outlined by Reference Zienkiewicz, Godbole, Gallagher, Gallagher, Oden, Taylor and ZienkiewiczZienkiewicz and Godbole (1974) for extrusion problems. This type of formulation, however, requires use of plate-bending type elements with C ₁ continuity which demands significant computational resources. Furthermore, it is difficult to satisfy stress-dependent sliding boundary conditions and to model compressible flow fields. In the approach proposed herein, velocities are first computed by using a primitive variable formulation, and then particle paths are determined through a stream-function formulation that makes use of the kinematic relationship

(1)

where ѱ is the stream function and ω is the vorticity of the flow field. The vorticity is calculated in the primitive variable formulation by using the relationship

(2)

where v ₁ is the velocity in the x ₁ direction and v ₂ is the velocity in the x ₂ direction. Equation (1) is obtained by substituting into Equation (2) the stream-function gradients for velocities v ₁ = ∂ѱ/∂x ₂ and v ₂ = ∂ѱ/∂x ₁. This relationship has been used previously to locate free surfaces in polymer-extrusion problems (Reference MitsoulisMitsoulis, unpublished).

Primitive Variable Formulation

Reference StolleStolle (unpublished) developed a finite-element model for studying steady-state planar glacier flow subject to gravity loading. It is assumed in this model that the influence of elastic strains is negligible, the flow field is incompressible, and flow parameters can be expressed in terms of non-Newtonian fluid rheology.

The primitive-variable approach adopted for this model is based on the principle of virtual velocities, which is used to write the equilibrium equations in a suitable integral form

(3)

where δv _i and δd _ij are virtual velocities and compatible virtual rates of deformation respectively, b_i are body forces, V is the volume, and T _i are surface tractions on the boundary r. The stresses σ_ij are related to the deformation rates d_ij and mean normal stress σ_m by

(4)

with

(5)

where µ is the deformation-rate and temperature-dependent viscosity given by

(6)

The symbol δ_ij is the Kronecker delta, σ_e = (3/2S _ij S _ij)^1/2 and = (2/3d _ij d _ij)^1/2 are Dorn’s definitions for equivalent stress and deformation rate, respectively, and S _ij are deviatoric stresses. Summation over repeated indices is assumed. A creep relationship between effective deformation rate and stress is required. A power-law form, generally used in glaciology, has been adopted for this paper:

(7)

where A is a temperature-dependent scalar function and n is a constant. At this point it should be noted that from a practical point of view no distinction need be made between strain-rate tensor generally used when describing material behaviour derived from ice-deformation tests and rate of deformation tensor which is associated with the principle of virtual velocities, provided that one is dealing with small strains.

A second virtual work-rate equation is necessary to enforce incompressibility

(7)

This equation states that the internal work rate due to a virtual mean normal stress δσ_m, acting on an incompressible flow field is zero.

The primitive variables are mean normal stress, σ_m, and velocities, v _i. A similar finite-element model has been developed by Reference Hooke, Hooke, Raymond, Hotchkiss and GustafsonHooke and others (1979), in which they used quadilateral finite elements. They observed longitudinal oscillations in predicted pressure which they attributed to their model. To avoid such oscillations, quadratic velocity and linear mean normal stress interpolations were used within each six-noded triangular element. Four-point integration was used to assemble the stiffness matrix and load vector. For details on the assembly of the finite-element equations, the reader is referred to Reference ZienkiewiczZienkiewicz (1977).

Stream-Function Formulation

The Galerkin weighted residual procedure was used to develop the stream-function finite-element equations. The weighted residual expression corresponding to Equation (1), subject to natural type and essential boundary conditions is given by:

(8)

where δѱ is the weighting function consistent with ѱ, v _s is the tangential boundary velocity, and the other terms the same as defined previously. Again, the reader is referred to Reference ZienkiewiczZienkiewicz (1977) or any other finite-element method text for details on the reduction of Equation (8) to a matrix form that is suitable for computer implementation.

The stream-function finite-element model was subjected to numerous tests to verify its accuracy. Originally, a three-noded triangular with linear interpolation for stream function was employed. However, the element performed poorly and was replaced with a six-noded triangular element with quadratic interpolation. Input for the stream-function finite-element program from the primitive variable simulations consisted of vorticity at each integration point and boundary velocities. The zero streamline was specified to conform to the bedrock topography.

To aid interpretation of results from the stream-function model, a computer program was developed to draw contours of constant stream function. For an incompressible flow field, these contours represent the paths that particles follow as the ice creeps. This program also calculated the resident times of ice particles at points along particle paths which could then be used to draw the isochrones.

Examples

Idealized double-slope ice mass

The primitive variable and stream-function finite-element models were used to study an idealized ice mass, the double-slope ice mass (Fig. 2 ). Flow parameters obtained by Hooke for the Barnes Ice Cap were used to simulate the creep behaviour of the double slope:

(9)

where σ _e is measured in bars and ϵ _e is in units of per annum. A unit weight of 8.952 kN/m³ was assumed. It should be noted that the constant in the flow law is consistent with units and definition of effective stress adopted.

Fig. 2. Double-slope ice mass idealization.

The boundary conditions for the creep simulations with the primitive variable formulation were a stress-free upper surface, a no-slip boundary along the bed, and only vertical movement along the divide. The boundary conditions for the stream-function simulations with the stream-function finite-element formulation were ѱ = 0 along the bed and divide, and velocities were specified along the upper surface of the ice mass. The results of the stream-function simulations are shown in Figure 3. It is clearly shown that the isochrones are approximately parallel to the ice-mass surface. A sensitivity study indicated that the use of a coarser grid, i.e. three elements, did not significantly affect the predictions of particle paths. In general, however, it may be said that the more elements there are, the better the representation of the actual continuum, and the more accurate the predictions for particle paths and isochrones. Of course, one must keep in mind the limitations of the mathematical models when attempting to use model predictions to assess actual field behaviour. That is, if the mathematical model used to describe the flow behaviour is poor, then any predictions provided by the model are likely poor. The variation of material parameters, used by the model, throughout the continuum must also be properly taken into account.

Fig. 3. Particle paths and isochrones for double-slope problem.

Barnes Ice Cap

The Barnes Ice Cap is a medium-sized ice cap located on Baffin Island, N.W.T., Canada. The geometry of the ice mass, along the Trilateration Net Flowline was taken from Reference Hooke, Hooke, Raymond, Hotchkiss and GustafsonHooke and others (1979); the finite-element grid used in this study is shown in Figure 4. Two types of ice exist in this ice mass at this section: white ice, located near the base, and blue ice which composes most of the ice cap. While the difference in unit weights of each ice was accounted for in the ice-dynamic analysis, it was assumed that both types of ice are isothermal and incompressible.

Fig. 4. Barnes Ice Cap idealization: (a) variation of A parameter; (b) grid.

The assumption of incompressibility for the entire problem is believed to be reasonable since the ice cap is nourished “by the formation of superimposed ice rather than snow-firn-ice metamorphism" (Reference ClassenClassen, 1977). Since this paper focuses on techniques of analysis, no attempt has been made at this point in time to account for actual temperature variations which influence the creep law. In order to establish reasonable ages for glacier ice, the velocities generated by the primitive variable finite-element model must be close to those measured in the field. To obtain a reasonable match between predicted and measured horizontal surface velocities, the A parameter of the creep law was adjusted along the glacier, while the power was held constant. Since more detailed information on velocities at depth were not available, no attempt was made to vary the A parameter with depth. At this point, it should be noted that the optimization procedure accounts for mean variations in temperature and fabric along the ice mass via the A parameter. The increase in A parameter toward the margin as shown in Figure 4 likely reflects increased thickness of softer white ice relative that of blue ice. The boundary conditions for the primitive variable finite-element program were a no-slip boundary at the base, a stress-free upper surface, and only vertical deformation at the ice-cap divide. For the stream-function finite-element simulation, ѱ = 0 was imposed on the bed and along the divide, and velocities were specified along the upper surface.

Figure 5 compares computed and measured surface velocities. While reasonable agreement exists between horizontal velocities, computed vertical surface velocities are consistently higher than the measured ones, except at one point. Reference Hooke, Hooke, Raymond, Hotchkiss and GustafsonHooke and others (1979) attributed these discrepancies to the influence of finite-element discretization on solution, errors associated with describing ice-cap geometry, and the effect of small transverse strain-rates. Grid error is responsible for some of the difference. However, the quasistatic equilibrium analysis completed in this study is expected to underestimate vertical surface velocities since the influence of fabric and vertical temperature variations, which would concentrate shear closer to bedrock thereby increasing uniformity of longitudinal strain-rates with depth, was not taken into account. Furthermore, the effect of transverse strain-rates is considered to be negligible for this section since the deviation of the principal strain-rates from the flow line is less than 3 (arc) as reported by Reference HoldsworthHoldsworth (1975). Thus, the authors believe that a major factor contributing to the lack of agreement between computed and measured velocities is the error associated with defining the ice-cap geometry.

Fig. 5. Surface velocities for the Barnes Ice Cap: (a) horizontal; (b) vertical.

The resulting particle paths and isochrones for the Barnes Ice Cap are shown in Figure 6. It can be clearly seen that the particle paths are significantly affected by changes in the bedrock topography and that particle paths in the upper part of the glacier are affected as well as those near the bedrock. It should be noted that this observation is consistent with that of Reference Reeh, Reeh, Johnsen, Dahl-Jensen, Langway, Langway, Oeschger and DansgaardReeh and others (1985), who recently showed that both isochrone and flow-line predictions are sensitive to bedrock topography.

Fig. 6. Particle paths and isochrones for the Barnes Ice Cap problem.

The location of the 10000 year isochrone is of particular interest. Reference HookeHooke (1976) has established that the top of the white ice at the base of the Barnes Ice Cap is of Pleistocene age, that is approximately 10000 years old. The 10000 year isochrone determined from the finite-element analyses approximately matches the top of the white ice, thereby confirming, by independent means, Hooke’s assertion.

Conclusions

An attempt has been made to define better the velocity field for purposes of refining particle-path and isochrone predictions by using the finite-element method. While simplifying assumptions regarding ice rheology and temperature distribution were assumed for the analyses presented in this paper, the finite-element method has the flexibility to accommodate further refinements in analyses if more detailed input information is available.

From this study two main conclusions can be reached:

• (1) The finite-element method procedure outlined in this paper has great potential for accurately determining the age of glacier ice, provided that the calculated and measured velocities almost match. A major drawback in the procedure lies in the manner in which the calculated velocities are fine-tuned to correspond to measured surface velocities. In the absence of field data regarding velocities at depth, further refinements to the procedure will not guarantee improvements in predictions.

• (2) Particle paths even near the surface of the Barnes Ice Cap seem to reflect changes in bedrock topography.