A framework for the prediction of solute spread in partially saturated heterogeneous porous formation is presented by combining a general Lagrangian formulation, relating the spatial moments of the solute body to the velocity field, with the stochastic theory of Yeh et al. (1985a,b) for steady unsaturated flow, relating the statistical moments of the velocity field to properties of the heterogeneous formation. The relevant formation properties are log saturated conductivity, log Ks, and the pore size distribution parameter, &agr;, which, in turn, are viewed as independent, random spatial functions. First-order approximations of the longitudinal and transverse components of the velocity covariance function are derived for unidirectional, vertical mean flow in partially saturated, heterogeneous porous formations of three-dimensional structure. Neglecting pore scale dispersion, and by using basic, kinematical relationships to relate the displacement to the velocity field, the effect of mean water saturation, and of statistical parameters of porous formation properties on longitudinal and transverse components of the time-dependent displacement covariance tensor X, is evaluated. Results of the analyses show that solute spread increases as water saturation decreases. In the case where the correlation scale of &agr; was small, as compared with that of log Ks, however, the longitudinal component of X is essentially robust to water saturation, suggesting that in this case, vandose zone transport might be amenable to an analysis similar to that for saturated groundwater transport. ¿ American Geophysical Union 1993