Statistical moments of solute plumes from small sources in variably saturated, heterogeneous porous media are analyzed by using a newly developed, efficient high-resolution Monte Carlo technique. In agreement with previous theoretical work, it is illustrated that the prediction of such solute plumes is associated with large uncertainties for dimensionless travel times, t', exceeding 40, particularly predictions of plumes in highly heterogeneous soils (&sgr;2y>2). Uncertainty about the travel path of the plume center contributes significantly to overall concentration uncertainty as flux fields become more variable. It is shown that the concentration coefficient of variation at the center of the plume initially increases but stagnates or decreases at later times. For highly heterogeneous soil flux conditions or for the common case of soils with strongly anisotropic conditions, analytical models underestimate transverse spreading of the mean concentration plume at any given time, while overestimating longitudinal spreading. At identical mean plume displacement distances, analytical models underestimate both transverse and longitudinal spreading and overestimate the variance of solute flux (breakthrough curve) by up to a factor 4. As an alternative to the statistical analysis of solute flux, we propose to analyze statistical properties of time associated with peak solute flux and with first exceedance of a given solute flux level. ¿ American Geophysical Union 1996