Particle tracing in ice sheets is required to deal with problems such as ice dating, oxygen isotope contents, or the distribution of any conservative property that is advected with the ice. The Lagrangian approach, in which a particles trajectory is constructed by numerical interpolation as it moves through an evolving ice sheet, is conceptually straightforward, but demanding in terms of its practical implementation. The main advantage of the method as compared to a Eulerian approach is that it is diffusion free, and that it immediately yields the trajectories of particles and the distribution of transported properties on isochronous surfaces. Its optimal implementation requires an accurate balance between computational overhead and desired accuracy. Because the Lagrangian method consists of the consecutive computation of a position of a set of particles within a given grid, the central problem is one of interpolation. Its accuracy depends on the interval and initial spatial distribution of particle launches, but generally leads to huge volumes of data to be stored and processed. We have implemented a Lagrangian tracer algorithm in a time-dependent thermomechanical model of the Antarctic ice sheet. The model has components describing ice-sheet and ice-shelf flow, isostatic adjustment adjustment of the lithosphere, and the derivation of past environmental boundary conditions. Numerical experiments are carried out for the last 4 glacial cycles forced by the Vostok temperature record. We will examine the time and spatial variations of two characteristics crucial for palaeoclimatic studies ice age (time of deposition) and the oxygen isotope content. These results are also compared to those obtained by solving the advection equation in the Eulerian way, and advantages and disadvantages of both approaches are discussed.