In this paper a stochastic model is developed to predict the behavior of partitioning tracers in a three-dimensional aquifer subject to spatial heterogeneities of nonaqueous phase liquid (NAPL) and Darcy flux. Stochastic Eulerian perturbation methods are applied to derive a system of six coupled partial-differential first-order moment equations from the traditional advection-dispersion-retardation equation. These equations describe the transient behavior of the concentration moments that explicitly accounts for the effects of the considered heterogeneities on prediction uncertainty and spatiotemporal correlation between tracer concentration, Darcy flux, and NAPL distribution. A finite difference technique is used to solve these equations numerically. Three different cases are presented to evaluate the impact of various statistical relationships between the random fields of natural log conductivity and natural log volumetric NAPL content. It is found that strong negative correlation between NAPL content and hydraulic conductivity leads to increased prediction uncertainty, whereas strong positive correlation between NAPL content and hydraulic conductivity leads to decreased prediction uncertainty over the uncorrelated case. For all cases, cross correlations between tracer concentration and Darcy flux or NAPL saturation are more extensive in the direction of mean flow than transverse to the mean flow. Furthermore, the magnitude of these correlations is approximately constant over the duration of the tracer experiment, implying that all concentration measurements taken at a particular location contain approximately equal information about Darcy flux or NAPL distribution within the experimental domain. The cross correlation between tracer concentration and Darcy flux is relatively insensitive to the degree of the conductivity-NAPL correlation. However, the correlation between tracer concentration and NAPL concentration is much larger for strongly (positively or negatively) correlated hydraulic conductivity and NAPL distributions. In a comparison paper <Zhang and Graham, this issue> these correlations, together with tracer measurements, are incorporated into a parameter estimation algorithm to produce site-specific conditional estimates of NAPL and Darcy flux distribution. ¿ 2001 American Geophysical Union