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...