A framework for the modeling of solute flux and related entities in partially saturated, heterogeneous porous formations was presented by combining a general Lagrangian formulation (Dagan et al., 1992), relating the travel time moments of the solute pulse 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. First-order approximations of the travel time covariance were derived for unidirectional, vertical mean flow in partially saturated, heterogeneous porous formations of three-dimensional structures. Hypothesizing lognormal distribution for the one- and two-particle travel time probability density functions, the effect of mean water saturation and statistical parameters of the porous formation properties on expected values, and variances of the solute discharge S and accumulated mass M passing through a given horizontal control plane located at an arbitrary vertical distance from the solute source, as a function of time were evaluated. Assessment of the uncertainty in the predictions of S and M and possible application of the results of this investigation to assessment of groundwater contamination hazards were discussed. ¿ American Geophysical Union 1993