The temperature can be solved for in a simple 1-D column if we
know the vertical (
) and horizontal (
)
velocity distributions, as well as the horizontal temperature
gradients (
), and any internally
generated heat terms (
). In an equation of this from,
the only neglected term is the horizontal heat transport due
to horizontal conduction, and this is negligible.
The heat equation is just the following:

Expressing the total time derivative (
) in terms of a
partial time derivative and the vertical and horizontal
advection terms we obtain the following expression:

In solving this equation the z direction terms will be handled implicitly, while the x direction terms (horizontal velocity and temperature gradients) will be included with the internal heat generation term. These gradients will need to be derived from some external 1-D or 2-D grid in the horizontal direction.
Having obtained a solution for
we can then obtain the
ice hardness as a function of z from an expression of this
form:

The strain rate is expressed in terms of a flow law, and is itself the vertical gradient of the horizontal velocity, given by the following expression:

The stress,
, will be zero at the top surface
and equal to the driving stress at the bed. Lacking any better
models we will assume a linear distribution through the
thickness, given by the following expression:

Hence we obtain the following expression:

which we can integrate from the bed (where the velocity is just equal to the sliding component of the velocity) to some depth, z, at which we wish to know the velocity. We will express the total horizontal velocity in terms of a parameter f, which is the fraction of the velocity due to sliding.

The sliding velocity can be expressed by some sliding relationship such as the one derived by Weertman [].

The horizontal velocity at any depth, z, will just be given by the following expression:

The integral can be evaluated numerically.
For a column-integrated model we will also need the column-averaged velocity, which is obtained by integrating the velocity from the bed to the surface and dividing by the thickness.

This also can be evaluated numerically.
For the temperature solution we also need the distribution of the vertical velocity in the z direction. This can be obtained from the incompressibility condition which requires that the sum of the principle strain rates be zero. Expressing the strain rates as velocity gradients we have the following expression:

or

In a flowband (1-D) model, the y term would be expressed in terms of the converging or diverging width of the flowband. In a map-plane (2-D) model, the horizontal velocity would have two components (x and y), each of which would be obtained from a treatment similar to the one above from which we obtained the x component of the velocity.
This can be solved for the gradient of the vertical velocity with respect to z and integrated from the bed to some depth z, obtaining the following expression:

Again, this integral can be evaluated numerically. We can probably assume that the vertical velocity at the bed is zero, in the absence of melting or freezing.
At this point we have both the vertical and horizontal velocities as a function of depth as required for the temperature solution. We can repeat the solution for this new distribution until the overall solution is self-consistent.