Eulerian models developed to simulate dispersion in fluid mechanics often consider the flux of the contaminant species to be proportional to the concentration gradient via a constant or time-dependent dispersion coefficient. These models are crude approximations for systems with velocity fluctuations evolving over a hierarchy of scales on the scale of observation. We say a system behaves in a Fickian fashion if the dispersion coefficient is constant, it is quasi-Fickian if the dispersion coefficient is time dependent, and it is convolution-Fickian if the flux is a convolution. The fractional flux in the sense of fractional derivatives is a special case of a convolution-Fickian flux. More general forms of the flux are possible, and in any case we call all fluxes anomalous if there is not a constant coefficient of proportionality between the flux and the gradient of concentration. In paper I of this two-part sequence we present a theory with statistical mechanical origins for simulating anomalous dispersion. Under appropriate limiting conditions the theory gives rise to Fickian, quasi-Fickian, convolution-Fickian, and fractional-Fickian fluxes. The primary result is a dispersive flux of integral type which in its most general form is not a convolution on time (it is non-Markovian however), but it is always a convolution in space. The concentration is represented by the inverse Fourier transform of the self-part of the intermediate scattering function. In paper II we present an experimental procedure that uses this theory to examine if and when the Fickian limit is reached in porous media homogeneous on the Darcy-scale but heterogeneous on the pore-scale. (C) 2001 American Institute of Physics.