A thermomechanical flow model has been developed in order to perform simultaneous calculations of the surface elevation and the distribution of stresses, velocities, and temperatures along two-dimensional flow lines on large grounded ice sheets. In order to solve the complex system of full thermomechanical coupled equations the following approximations have been made: A coordinate scaling giving a lead order reduced model and an iteration procedure, which decouples the energy balance equation from the rest of the equations. Glen's flow law is used as the constitutive relationship between stresses and deformation in the model, and it is found that the longitudinal deviatoric stresses have a significant role in the lead order equations. At the ice divide the solution shows that the velocity, stress, and temperature distributions change rapidly. Here the basal temperature increases creating a ''hot spot.'' The surface strain rates increase 50% and the horizontal velocity profile is concave near the base, with an inflection point. Away from the ice divide, basal temperature increases, and in the example presented here the basal temperature reaches the pressure melting point after 2/3 of the lateral extent and basal sliding begins to influence the flow. In the ablation region it is seen that the heat produced by internal deformation is of the same order of magnitude as the geothermal heat flux from the bedrock. ¿ American Geophysical Union 1989