Introduction

Previous investigators have solved a variety of glacier-flow problems. However, most of the tractable problems have had either simple geometries and boundary conditions or simplified constitutive equations for ice. The object of the research described in this paper is to determine whether direct methods from the calculus of variations can provide convenient approximate solutions of relatively complex glacier-flow problems.

Solutions of two sample problems are presented. An analysis of the closure of an intergranular vein in ice at the melting point is based upon a variational principle for velocities. An analysis of the flow of ice in a cylindrical channel of parabolic cross-section is based upon a variational principle for stresses.

Closure of an intergranular vein

A vein forms along each intersection of three crystals in ice at the melting point and the network of veins in a glacier provides a drainage system for water. If the overburden pressure in the ice is balanced by the veinous water pressure and the surface tensions of the curved ice-water interfaces, the vein will be in equilibrium. If the veinous water pressure is somehow reduced, deformation of the ice will tend to close the vein. A variational principle for velocities presented by Reference JohnsonJohnson (1960) can be used to determine, the closure rate.

The equilibrium shape of an intergranular vein has been described by Reference Nye and FrankNye and Frank (1973). If the ice is assumed to be isotropic and if the stress at some distance from the vein is assumed to be isotropic, the cross-section has the three lines of symmetry shown in Figure 1. Walford measured the angle subtended at each vein corner and reported a value of about 33⁰ to Reference Nye and MaeNye and Mae (1972).

Fig. 1. Shape of vein and finite-element model.

The vein is assumed to be cylindrical, and the variation of veinous water pressure along the vein is neglected, so the solution does not depend upon a coordinate parallel to the axis of the vein. The cross-section has three lines of symmetry, so the behavior of the ice between adjacent lines of symmetry typifies the closure of a vein.

The computational techniques developed for the vein-closure problem are detailed in a report by Reference OakbergOakberg (1978), which is available upon request and which has also been deposited with World Data Center A for Glaciology, Boulder, Colorado, U.S.A. Briefly, the Ritz technique is used to minimize the appropriate form of Johnson’s variational principle for velocities. A linear combination of coordinate functions is constructed, and the functional of Johnson’s variational principle is minimized with respect to the coefficients of the individual coordinate functions.

Two types of coordinate function are employed. One coordinate function is the velocity distribution that would result from the closure of an equivalent circular vein. This coordinate function reflects the radial symmetry of the velocity distribution that is expected at some distance from the vein. The other coordinate functions are constructed from the finite-element model shown in Figure 1. The coordinate functions consist of velocity distributions that are linear over each finite element. There are 152 nodal velocity components, and 19 of the nodal velocity components are specified. There are 123 elements in the model, and the area of each element is required to remain constant during deformation. Therefore, 123 constraint equations, in addition to the 19 specified velocity components, reduce the number of independent nodal velocity components to ten. These ten independent velocity components serve as variables in the minimization process, along with the radius of the equivalent circular vein.

The two types of coordinate function play different roles. The radially symmetric coordinate function establishes the solution at some distance from the vein. The finite-element coordinate functions adjust the radially symmetric solution in the vicinity of the vein.

Glen’ s formula as generalized by Reference NyeNye (1957) provides a constitutive equation for ice, and a value of 3 is used for the exponent. The velocity distribution on the ice-water interface is shown in Figure 2. The velocity distribution at some distance from the vein displays the expected radial symmetry.

Fig. 2. Non-dimensional nodal velocities on the ice-water interface, multiplied by 10000.

Reference NyeNye (1976) has employed a vein of circular cross-section that has the same closure rate as an actual vein. The radius of such an equivalent vein is 0.508 times the radius that circumscribes the actual vein.

As a check on the analysis, the stress distribution can be evaluated and compared with physical intuition. A least-squares technique is used to evaluate the nodal rate of deformation components, and the stress components are evaluated from the components of the rate of deformation. Isoclinics of the stress tensor, obtained by interpolation between nodal orientations, are shown in Figure 3. There are three isotropic points, one at the vein corner, one at infinity, and one along the ice-water interface.

Fig. 3. Isoclinics of the stress tensor.

The stress at each vein corner is isotropic. This is perhaps surprising, since the vein corner is re-entrant. An analysis by Reference OakbergOakberg (1978) shows that a singular solution would require a concentrated force at each vein corner, which could be provided by surface tensions that are not in equilibrium. However, such a solution would require a temperature distribution that is singular at each vein corner, and this rules out the singular solution on physical grounds.

The overall equilibrium of the model can be checked. First, the resultants of the tractions on the boundaries of the model are evaluated. The system of forces is very nearly concurrent but the magnitude of the resultant of the system is about 10% of the magnitude of an individual force. The lack of force equilibrium can be traced to a probable underestimation of the stress near the vein corner. The stress gradients are large near the vein corner, and the finite-element mesh is apparently too coarse in this region.

Flow of ice in a cylindrical channel of parabolic cross-section

If a glacier is not frozen to its bed, its flow is due to slip between the basal ice and the bed and to deformation of the ice. Most treatments of the problem have considered the deformation and the basal slip separately, but a variational principle for stresses presented by Reference JohnsonJohnson (1960) can be used to investigate the interaction between the two modes of flow.

Glen’s formula as generalized by Reference NyeNye (1957) provides the constitutive equation for ice, and a value of 3 is used for the exponent. The velocity–traction relationship developed by Reference WeertmanWeertman (1957) models the basal slip.

The computational details for the ice-flow problem are very similar to those for the vein-closure problem. The main difference is that the coordinate functions represent the state of stress within the ice instead of the velocity distribution.

Two types of coordinate function are employed. The first satisfies the differential equation of motion, the specified traction boundary conditions, and a criterion for stress gradients near the center-line of a glacier presented by Reference RaymondRaymond (1974). The other coordinate functions are generated from the finite-element model shown in Figure 4. The coordinate functions consist of shear-stress distributions that are linear over each finite element. There are 26 nodes in the model, so there are 52 nodal shear-stress components. Six traction components on the glacier surface and six traction components on the plane of symmetry are zero. There are 34 elements and the shear-stress components over each finite element are required to satisfy the homogeneous differential equation of equilibrium. These 34 equations, in addition to the 12 specified traction components, reduce the number of independent nodal shear-stress components to six. These six independent shear stress components serve as variables in the minimization process.

Fig. 4. Finite-element model and non-dimensional nodal velocities, multiplied by 1000.

The two types of coordinate functions play different roles in the solution. The basic coordinate function represents the stress distribution at some distance from the glacier bed, while the finite-element coordinate functions adjust this distribution near the glacier bed. The solution is presented in Figure 5. The contours of maximum shear stress are obtained by interpolation between nodal values.

Fig. 5. Contours of maximum shear stress.

The velocity distribution is obtained in the following manner. Nodal components of rate of deformation are obtained from the shear-stress components. The components of rate of deformation are integrated over each finite element to provide values for the relative nodal velocities. In addition, the basal slip velocities are evaluated from the velocity–traction relationship at the glacier bed. The resulting equations are overdetermined, and the nodal velocity distribution that is best in the least-squares sense is calculated. The nodal velocity distribution is presented in Figure 4. This particular velocity distribution has a relatively large contribution from basal slip.

As a check, the velocity distribution for a glacier without basal slip can be obtained and compared with a solution presented by Reference NyeNye (1965). The solutions agree to the three significant figures quoted by Nye.

Assessment of the methods

The generation of coordinate functions from finite-element models is not an essential feature of the methods. However, the author has not been able to discover other ways of obtaining coordinate functions that are suitable for complex problems. In the sample problems, coordinate functions obtained from finite-element models are combined with coordinate functions that are exact solutions of related problems. The finite-element coordinate functions make localized adjustments to the basic solutions.

Solutions based upon the variational principle for stresses are more satisfactory than solutions based upon the variational principle for velocities. In the vein-closure problem, the stress field is obtained from the velocity distribution. The process is essentially one of differentiation, which increases the relative errors of the distribution. A great deal of smoothing and careful interpretation is required to produce a realistic stress distribution. In the glacier-flow problem, the velocity distribution is obtained from the shear-stress distribution. The process is essentially one of integration, which decreases the relative errors of the distribution.

It is not always possible to construct coordinate functions from finite-element models: a velocity distribution can be overconstrained by the incompressibility requirements for the finite elements and the specified velocity boundary conditions; a stress distribution can be overconstrained by the equilibrium requirements for the finite elements and the specified traction boundary conditions. For example, it would not be possible to use the variational principle for stresses in the vein-closure problem, since the stress distribution would be overconstrained.

The functional and the coordinate functions that are employed depend upon the particular problem being solved. Therefore, it is unlikely that programs having general application could be developed. A great deal of programming effort is required for each problem.