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.
