A review is given of the theory of polythermal ice sheets, ie, ice masses of which the ice has submelting temperatures in certain subdomains and is at the pressure melting point in other subdomains. Cold ice is treated as a non-linearly viscous heat conducting fluid, temperate ice as a mixture of ice at the melting point and melt-water diffusing through the ice matrix. Cold and temperate ice are separated by a non-material singular Stefan-type surface. We repeat and partly amend the complicated field equations and boundary conditions as derived in the literature. These equations are subjected to a scale analysis that makes the creeping flow conditions and the shallow geometries of land-based ice sheets explicit. The small aspect ratio ε — typical depth to horizontal distance over which the geometry and/or stresses change appreciably — suggests a perturbation approach for a possible analytical or numerical solution which has been pursued to include second-order terms O (ε2 ). The lowest-order O (ε0 ) model equations, known as the shallow-ice approximation (SIA), are asymptotically valid in the entire ice sheet domain except a small marginal zone provided the topographic variations are shallow, ie, possess wave height to wavelength ratios that are O (ε1 ) and the constitutive relation for the stress deviator exhibits finite viscosity at zero effective shear stress (square root of second stress deviator invariant). We critically review earlier procedures and put them into the proper perspectives with regard to the original expansion procedures. We extend the zeroth-order theory to first and second order but present only those equations and deductions from them which lead to improved physical insight. In particular we derive stress formulas which show how the stresses depend on i) depth and surface slope, ii) surface topography and iii) stress deviator components, more complete than, and going beyond, known formulas of the literature. Finally we discuss numerically computed second-order stresses for the present state of the Greenland ice sheet. It turns out that they are typically three orders of magnitude smaller than the corresponding zeroth-order quantities, and that they are mainly determined by contributions due to zeroth-order stress deviators, rather than by topography effects. Their relative importance is largest close to the ice surface for the second-order pressure, and in the vicinity of ice domes for the horizontal, bed-parallel shear stresses. There are 229 references.