This paper presents a new Eulerian-Lagrangian method for modeling flow and transport of passive solutes in heterogeneous porous formations. The physically plausible random velocity fields are generated by a geostatistical based model. The ability of the method to correctly reproduce the velocity field statistics and to satisfy mass balance is tested and demonstrated for the case of two-dimensional flow. The method is found particularly useful in dealing with large fields where numerical, fixed grid methods such as finite difference or finite elements become very demanding from a computational point of view. The transport problem is solved in a Lagrangian framework by the particle-tracking approach which uses a suitable grid refinement to correctly handle the velocity fluctuation at any a priori defined resolution. The new method is employed to compute the probability distribution functions of the concentration and travel times and to investigate the limitations of existing methods for predicting concentrations. ¿ American Geophysical Union 1994