Abstract A method is described for simulating three-dimensional (3D) flow fields for use in transport modelling of contaminant migration in statistically anisotropic heterogeneous aquifers. This method differs from previous approaches in that the flow field is generated directly, reducing the computational requirements considerably by bypassing the interim steps of the more traditional techniques. The approach discussed here provides a more computationally efficient and verifiable means of modelling transport in a 3D formation. The method will enable the generation of concentration and travel time probability distribution functions (PDFs) as well as 3D images of the contaminant plume. INTRODUCTION Models which can simulate the migration of reactive solutes in a three dimensional (3D) flow field within anisotropic heterogeneous aquifers are necessary to more accurately estimate the movement of contaminants in the environment and to better understand the underlying processes which govern contaminant transport. Estimates limited to only the mean and variance of the solute concentrations are often insufficient since contaminant concentration distributions are typically non-Gaussian and not fully characterized by their first two moments (Smith & Schwartz, 1980; Bellin et al., 1994; Rubin et al., 1994). Until recently, the complete modelling of contaminant transport in heterogeneous porous media has been limited to analysis in two dimensions. Though traditional approaches such as finite element can be applied in three dimensions, the computational requirements are excessive. The method presented here is an alternative to these approaches and requires less computational effort without loss of accuracy. METHODOLOGY This method is an Eulerian-Lagrangian approach and is applied in a model utilizing a stochastic framework to predict the spatial moments and probability distribution functions (PDFs) of the solute concentration and travel time. The fluid velocity is