The flow of ice sheets and glaciers dissipates significant amounts of heat, which can result in the formation of ‘temperate ice’, a binary mixture of ice and small amounts of melt water that exists at the melting point. Many ice masses are polythermal, in the sense that they contain cold ice, below the melting point, as well as temperate ice. Temperature and melt water (or moisture) content conversely affect the flow of these ice masses through their effect on ice viscosity and sliding behaviour. Ice flow models therefore require a component that can solve for temperature and moisture content, and determine the free boundary between the cold and temperate subdomains. We present such a model, based on the theory of compacting partial melts. By contrast with other models, we describe gravity- and pressure-gradient-driven drainage of moisture, while maintaining a divergence-free ice flow at leading order. We also derive the relevant boundary conditions at the free cold–temperate boundary, and find that the boundary behaves differently depending on whether ice enters or exits the temperate region. The paper also describes a number of test cases used to compare with a numerical solution, and investigates asymptotic solutions applicable to the limit of small compaction pressure gradients in the temperate ice regions. A simplified enthalpy-gradient model is finally proposed, which captures most of the behaviour of the full model in this limit.