A conservative, semi-Lagrangian numerical method, domain of influence search for convective unconditional stability (DISCUS), for computing solute transport in fluvial systems with transient storage is presented. The model includes an explicit, shape-preserving cubic interpolation update for advection and fully implicit temporal treatments of dispersion and transient storage. Numerical results with this new method are compared against exact solutions in a uniform advective-dispersive flow with transient storage. Results from a conventional Eulerian numerical model are also shown. Both models are inherently stable and mass-conserving. However, the semi-Lagrangian method is shown to be accurate and free of grid-scale oscillations over a wide range of Courant and grid Peclet numbers in contrast to the Eulerian approach whose accuracy degrades at high Courant numbers and/or large-grid Peclet numbers. Some quantitative aspects of the comparative model efficiencies are discussed, with the semi-Lagrangian model being particularly attractive when large time steps are used. In contrast to earlier semi-Lagrangian approaches, DISCUS employs the use of a cumulative solute mass variable in the advective update, and the issue of its interpolation is examined in the paper. Linear interpolation is shown to guarantee bounded solutions but can be very inaccurate, with numerical errors manifesting themselves by introducing artificially enhanced longitudinal dispersion. Higher-order interpolation (cubic) increases accuracy but at the expense of introducing spurious maxima and minima which manifest themselves as small grid-scale oscillations. By implementing a flux-limiting technique, based on convexity preservation arguments, completely bounded solutions can be achieved in both spatial and temporal domains with the higher-order interpolation. ¿ 2001 American Geophysical Union