Accurate dating in ice sheets is required for a correct interpretation of palaeoclimatic records and for incorporation of material characteristics in the flow law which depend on ice age. In this paper we make a detailed comparison between a Lagrangian and Eulerian approach to the ice advection problem in numerical ice sheet models. This comparison is first performed for a schematic two-dimensional ice sheet of Nye-Vialov type with a prescribed stationary velocity field. Several cases are examined which incorporate basal melting, basal sliding, and an undulating bed. A further comparison is made with an analytical solution for the ice divide. Both methods are also applied in a thermomechanical model of the Antarctic ice sheet for present-day climate conditions. Our main conclusion is that for similar discretisation parameters the Lagrangian method produces less error than an Eulerian approach using a second-order upwinding finite-difference scheme, though the difference is small (<1%) for the largest part of the model domain. However, problems are introduced in the Lagrangian approach due to the dispersion of tracers necessitating the use of interpolation procedures that are a source of additional error. It is also shown that a cubic-spline approximation of Lagrangian trajectories improves accuracy, but such a method is computationally prohibitive in large-scale ice-sheet models.