Velocity covariances for saturated heterogeneous formations with anisotropic structures were derived for mean flow parallel and perpendicular to the formation bedding and were employed for modeling transport in these formations by the use of a Lagrangian formulation. Relatively simple analytical expressions, each requiring only a single numerical integration, were obtained for the longitudinal and transverse components of the velocity covariance, the displacement covariance, and the macrodispersion coefficient tensors, all as functions of the anisotropy of the formation. Closed-form analytical expressions were derived for the longitudinal and transverse components of the velocity variance and for the asymptotic values of the transverse components of the displacement covariance. The applicability of this approach for modeling flow and transport in the vadose zone is tested in the companion to this paper by the use of the stochastic theory of Yeh et al. (1985a, b) for unsaturated flow. ¿ American Geophysical Union 1995