In this work, three finite volume discretisations of the two-dimensional diffusion equation, governing groundwater flow and for use with structured quadrilateral meshes, have been developed. The three methods rely on a cell-centred finite volume approach, but show distinct differences in the choice of: gradient approximation, head interpolations and control volume. They have been implemented into MODFLOW to enhance its capability in simulating accurately natural aquifer boundaries. A time implicit formulation has been used in each model. Five test examples have been undertaken to compare the performance of the newly developed methods against MODFLOW predictions and analytical results. The accuracy of the results obtained was found to depend on the spatial and temporal discretisations. One of the three developed methods proved its robustness, with regard to mesh non-orthogonality and skewness, and has been tested with a field case study. A discussion about the performance of the new developed finite volume model has been included and the model has been shown to perform well in comparison with MODFLOW.