A new three-dimensional thermomechanically coupled ice sheet model is developed. Contrary to the majority of three-dimensional ice sheet models (shallow ice approximation), higher-order stress gradients (longitudinal and transverse) are accounted for in the force balance equations. The horizontal velocity field is determined from the force balance equations in their “derivative form” (elliptical equation) using the method presented by Pattyn [2000 , 2002a]. The model is solved on a regular grid in the horizontal and an irregular grid in the vertical and is numerically stable. Basic experiments include the European Ice Sheet Modeling Initiative (EISMINT) benchmarks for large-scale ice sheet models and a comparison with the Saito-Blatter ice sheet model including higher-order stress gradients [ Saito et al., 2003 ]. Detailed calculations of ice flow over three-dimensional bedrock perturbations showed the validity of the higher-order solution. The model is capable of simulating the evolution of an ice stream within the ice sheet and shows important aspects of observed ice stream features, such as the surface flattening and the importance of side drag. The simulation of the ice flow over a subglacial lake results in a flattening of the surface, a local velocity increase over the lake, and a deviation of the ice flow from the main flow direction, features which are also observed at Lake Vostok, Antarctica.