Monte Carlo studies of flow and transport in two-dimensional synthetic conductivity fields are employed to evaluate first-order flow and Eulerian transport theories. Hydraulic conductivity is assumed to obey fractional Brownian motion (fBm) statistics with infinite integral scale or to have an exponential covariance structure with finite integral scale. The flow problem is solved via a block-centered finite difference scheme, and a random walk approach is employed to solve the transport equation for a conservative tracer. The model is tested for mass conservation and convergence of computed statistics and found to yield accurate results. It is then used to address several issues in the context of flow and transport. The validity of the first-order relation between the fluctuating velocity covariance and the fluctuating log conductivity is examined. The simulations show that for exponential covariance, this approximation is justified in the mean flow direction for log conductivity variance, &sgr;f2, of the order of unity. However, as &sgr;f2 increases, the relation for the transverse velocity component deviates from the fully nonlinear Monte Carlo results. Eulerian transport models neglect triplet correlation functions that appear in the nonlocal macroscopic flux. The relative importance of the triplet correlation term for conservative chemicals is examined. This term appears to be small relative to the convolution flux term in mildly heterogeneous media. As &sgr;f2 increases or the integral scale grows, the triplet correlation becomes significant. In purely convective transport the triplet correlation term is significant if the heterogeneity is evolving. The exact nonlocal macroscale flux for the purely convective case significantly differs from that of the convective-dispersive transport. This is in agreement with recent theoretical analysis and numerical studies, and it suggests that neglecting local-scale dispersion may lead to large errors. Localization errors in the flux term are evaluated using Monte Carlo simulations. The nonlocal in time model significantly differs from the fully nonlocal model. For small variance and integral scale there is a slight difference between the fully localized flux and the fully nonlocal convolution flux. This is also in agreement with recent theories that suggest that moments through the second for the two models are nearly identical for conservative tracers. The fully localized model does not perform well in the purely convective cases. ¿ 1998 American Geophysical Union