Numerical simulation of solute transport in heterogeneous porous media is greatly complicated by the large velocity and concentration gradients induced by spatial variations in hydraulic conductivity. Eulerian-Lagrangian methods for solving the transport equation can give accurate solutions to heterogeneous problems if their interpolation algorithms are properly selected. This paper compares the performance of four Eulerian-Lagrangian solvers that rely on linear, quadratic, cubic spline, and taut spline interpolators. In each case a tensor product decomposition is used to reduce the general n-dimensional interpolation problem to a sequence of n one-dimensional problems. Comparisons of a set of test problems indicate that the linear and taut spline interpolators are dispersive while the quadratic and cubic spline interpolators are oscillatory. The cubic and taut spline interpolators give consistently better accuracy than the more conventional linear and quadratic alternatives. Simulation experiments in two- and three-dimensional heterogeneous media indicate that the taut spline interpolator, which is applied here for the first time to a solute transport problem, is able to yield accurate essentially nonoscillatory solutions for high grid Peclet numbers. The cubic spline interpolator requires significantly less computational effort to achieve performance comparable to the other methods. ¿ 1999 American Geophysical Union