The numerical simulation of groundwater flow through heterogeneous porous media is discussed. The focus is on the performance of a parallel multigrid preconditioner for accelerating convergence of conjugate gradients, which is used to compute the pressure head. The numerical investigation considers the effects of boundary conditions, coarse grid solver strategy, increasing the grid resolution, enlarging the domain, and varying the geostatistical parameters used to define the subsurface realization. Scalability is also examined. The results were obtained using the ParFlow groundwater flow simulator on the CRAY T3D massively parallel computer.