Due to their low cost and high expansion capabilities, distributed memory parallel computers, consisting of relatively large arrays of processor elements with local memory, offer great potential for solving large scale three-dimensional subsurface flow problems. We consider here two possible solution methods for such detailed modelling using the finite element method. One is based on the commonly used diagonally scaled conjugate gradient method. The other uses an element-by-element preconditioner in the conjugate gradient method which exploits the natural parallelism in the finite element method. Results using both methods are documented for saturated and variably saturated flow simulations carried out on a single Intel i860 processor and an array of Inmos transputer/Intel i860 processors at Lancaster.