Modelling the theoretical response of several important geophysical systems involves the solution of Poisson's equation with homogeneous Neumann boundary conditions (i.e. a zero normal gradient) imposed over either open or closed surfaces.A simple integral equation solution to this problem is derived from first principles. It is applicable to both types of surface and in this respect represents an improvement on existing integral equation techniques. However, the present surface integral equation displays a strong singularity of order 1/R3 which requires an appropriate interpretation for its implementation.A comparison of some numerical results with analytical data taken from the literature demonstrates that the proposed integral equation technique is suitably robust, accurate and efficient for practical application in geophysical interpretation.