The flow of large natural ice masses under gravity is described by the mass, momentum, and energy balances of an incompressible, homogeneous, heat conducting, non-linearly viscous fluid in which the shear response includes a strongly temperature-dependent rate factor. Dimensionless analysis and co-ordinate stretching reflecting the long aspect ratio show that series expansions in a small parameter which determines the surface slope magnitude are uniformly valid even when temperature variation induces a strongly non-uniform mechanical response. The normalised energy balance shows that both horizontal and vertical advection are significant in thin and thick grounded sheets and in floating shelves, and that viscous dissipation can be significant in basal regions of a grounded sheet, and hence there is strong thermomechanical coupling. Moreover, though a thermal basal boundary layer may arise in a thick sheet, it would only give rise to significantly enhanced temperature and strain-rate gradients in extreme circumstances. The leading order relations for steady plane flow of a grounded sheet are reduced to a parabolic system for the temperature and two velocity components, which incorporates the unknown surface slope in coefficients and boundary conditions. This provides a useful starting point for numerical solution of the thermomechanically coupled problem. A fixed domain mapping is presented as an attractive alternative formulation when the bed topography is close to planar.