Particle tracing in ice sheets is required to deal with problems such as ice dating, oxygen isotope contents, or the distribution of any conservative characteristic 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 advantages of the method as compared to a Eulerian approach are 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. 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. Tracers are launched at the surface every 100 years on a grid with 20 km resolution. Their current positions are calculated using piece-wise linear interpolation. Two methods of numerical integration are examined the Eulerian one (not to be confused with the general Eulerian approach referred above) and the Runge-Kutta one at different time steps for calculating positions of tracers. The poster will display results for ice age (date of deposition) and isotopic composition. These results are also compared to those obtained by solving the advection equation in the Eulerian way.