[1] A mathematical model for polythermal glaciers and ice sheets is presented. The enthalpy balance equation is solved in cold and temperate ice together using an enthalpy gradient method. To obtain a relationship between enthalpy, temperature, and water content, we apply a brine pocket parameterization scheme known from sea ice modeling. The proposed enthalpy formulation offers two advantages: (1) the discontinuity at the cold-temperate transition surface is avoided, and (2) no treatment of the transition as an internal free boundary is required. Fourier's law and Fick-type diffusion are assumed for sensible heat flux in cold ice and latent heat flux in temperate ice, respectively. The method is tested on Storglaciären, northern Sweden. Numerical simulations are carried out with a commercial finite element code. A sensitivity study reveals a wide range of applicability and defines the limits of the method. Realistic temperature and moisture fields are obtained over a large range of parameters.Citation: Aschwanden, A., and H. Blatter (2009), Mathematical modeling and numerical simulation of polythermal glaciers,