A two-dimensional model is used to describe contaminant transport in the presence of colloids for heterogeneous porous media. The model accounts for both spatially varying conductivity (physical heterogeneity) and the spatially varying distribution coefficient and colloid attachment coefficient (chemical heterogeneity). The model is tested against experimental data, and the results are favorable. The model is implemented in a stochastic Monte Carlo fashion, and both absolute and relative dispersion frameworks are used for the analysis. One finding in this study indicates that the presence of colloids reduces variability in mass arrival times to a downstream control plane. The study also indicates that the effect of geochemical heterogeneity is important only if it is correlated to physical heterogeneity. Spatial variability of the contaminant distribution coefficient or colloid attachment coefficient (i.e., geochemical heterogeneity) is not important in determining the mean plume shape and concentration values when it is not correlated to the hydraulic conductivity variability. Moreover, the comparison between absolute and relative dispersion shows that including the plume meandering in the overall ensemble (absolute mean) would reduce the peak concentration 20--30%. It is also found that the positive correlation structure of geochemical parameters and hydraulic conductivity reduces the difference between absolute and relative dispersion results.