In the case of a non-isothermal glacier it is necessary to integrate the equations of dynamics together with the equation of heat conduction, heat transfer, and heat generation because of the interdependence (1) of strain-rate of ice on its temperature, and (2) of ice temperature on the rate of heat transfer by moving ice and on the intensity of heat generation in its strain. In view of the complexity of the whole system of equations, simplified mathematical models have been constructed for dynamically different glaciers. The present model concerns land glaciers with thicknesses much less than their horizontal dimensions and radii of curvature of large bottom irregularities, so that the method of a thin boundary layer may be used. The principal assumption is the validity of averaging over a distance of the order of magnitude of ice thickness.

Two component shear stresses parallel to the bottom in glaciers of this type considerably exceed the normal stresses and the third shear stress, so the dynamics are described by a statically determined system of equations. For the general case, expressions for the stresses have been obtained in dimensionless affine orthogonal curvilinear coordinates, parallel and normal to the glacier bottom, and taking into account the geometry of the lower and upper surfaces. The statically undetermined problem for ice divides is solved using the equations of continuity and rheology, so the result for stresses depends considerably on temperature distribution. In the case of a flat bottom the dynamics of an ice divide is determined by the curvature of the upper surface.

The calculation of the interrelating velocity and temperature distributions is made by means of the iteration of solutions (1) for the components of velocity from the stress expressions using the rheological equations (a power law or the more precise hyberbolic one) with the assigned temperature distribution, and (2) for the temperature with the assigned velocity distribution. The temperature distribution in the coordinate system used is determined by a parabolic equation with a small parameter at the principal derivative. Its solution is reduced to the solution of a system of recurrent non-uniform differential equations of the first order by means of a series expansion of the small parameter: the right part for the largest term of the expansion contains a function of the heat sources, and for the other terms it contains the second derivative along the vertical coordinate from the previous expansion term.

Thus advection makes the main contribution to the heat transfer, and temperature in a glacier is distributed along the particle paths, changing simultaneously under the influence of heat generation. A relatively thin conducting boundary layer adjoins the upper and lower surfaces of a glacier, playing the role of a temperature damper in the ablation area. The equation of heat conduction (at the free surface) or of heat conduction and heat transfer (at the bottom) with the boundary conditions, and with the condition of the connection with the solution of the problem for the internal temperature distribution, is being solved for the boundary layer because of its small thickness. Beyond the limits of the boundary layer, heat conduction makes a small change in the temperature distribution, which can be calculated with any degree of accuracy.