A discussion of the equations governing the flow of ice must begin with an evaluation of the conservation of mass. Mass conservation requires that the net flux into a region be balanced by changes in the ice thickness of that region. There are three factors involved in the mass balance of an ice sheet: 1) the rate of accumulation or ablation; 2) the change in ice surface elevation with respect to time; 3) the divergence of ice flux density.
Balancing fluxes into and out of a differential area with the rate of thickening in that area yields a continuity equation of the form
The accumulation rate, , represents the net total of both the snow fall rate and the melting rate at a point of the domain.
To complete this equation one must relate the flux, , to the ice surface elevation, h, through a constitutive equation. For non-Newtonian fluids such as ice, this involves a nonlinear relationship depending on the ice surface elevation, the thickness, and the surface slope. We can, however, linearize this by absorbing all dependencies on h and all nonlinearities in the surface slope into a spatially nonuniform proportionality constant, . This constitutive equation relating ice flux to surface slope can be expressed by
The material coefficient is dependent upon several material properties of ice. The physical basis for this relationship and the nature of the dependency of on the material properties of ice will be made clear with further evaluation of the ice velocity, U. The bed topography relates ice thickness, H, to ice surface elevation, h.
In glacier flow the ice behaves like a stack of pages, each page above moving relative to the one below it, with the bottom-most page affixed to the bed. The strain rate (change in length per unit length per unit time) is related to the driving stress by a non-linear flow law of the form [,,]
where the derives from this being a shear deformation. Neglecting longitudinal stresses, S at the bottom of the ice column is simply the basal traction, given by . Since the stress must be zero at the free top surface of the ice column, a simple expression for S as a function of depth is
where the coordinate system is such that z is zero at the bed and H at the surface. Integrating equation () from the bed where velocity is zero to some depth z one obtains the velocity as a function of depth .
Integrating equation () from the bed to the surface and dividing by the thickness yields the column-averaged ice velocity necessary for the flux calculation of the FEM.
Under certain conditions ice can slide over its bed. This process is not well understood, but evidence indicates that ice can move as a block with a nonzero velocity at the bed and with little or no internal shear. Sliding appears to depend on the concentration of debris in the ice at the ice-bed interface, the ice grain size, the presence of air bubbles in the ice, the permeability of the bed, the degree to which the the bed can be deformed, the amount of water at the ice-bed interface, the roughness of the bed, the subglacial water pressure, and cavitation .
The sliding law used here is a general relationship for beds at the melting point developed by Weertman [,]. It depicts movement over a rough bed which occurs as a result of melting on the high-pressure upstream side of an obstacle followed by subsequent refreezing on the low-pressure downstream side. The sliding law constant, B, and the sliding law exponent, m, can be used to absorb the effects of the various conditions influencing sliding. Weertman's theory yields the following expression for a sliding velocity that is uniform with depth.
A primary reason for incorporating the sliding law into an ice sheet model is to facilitate the study of ice streams. Ice streams are not well understood, but are believed to be caused by some loss of traction at the bed which allows sliding to occur.
The column-averaged ice velocity for a combination of these two modes of motion is
where f is the percentage of velocity due to sliding. For existing ice sheets f can be used as a parameter to `tune' the model. For reconstructing ice sheets from the past and for extrapolating future behavior f can be estimated from various erosional features found in the glacial geologic record .
The form of is obtained by combining equations (), (), and () and substituting into equation ().
This shows the dependence of the material coefficient on ice density, gravity, flow law constant and sliding law constant (both of which are dependent on ice temperature, ice crystal size, impurity content, etc.), ice thickness, and surface gradient.
Substituting the constitutive law expression for from equation () into continuity equation (), yields the following differential equation in h.
The presence of in the continuity equation incorporates the physics of the flow and sliding laws into the problem, since its form depends on the form of the flow and sliding law. Different treatments of the flow and sliding process change the form of but do not affect the method with which the problem is solved. An entirely different sliding relationship can be substituted without materially affecting the ensuing discussion leading to a solution to this equation.
With the continuity equation determined, all that remains for a complete description of the problem is a discussion of boundary conditions. The FEM allows for both essential and natural boundary conditions. Natural boundary conditions encompass situations where the flux is specified along a portion of the boundary, while essential boundary conditions require the specification of the ice height along a portion of the boundary.