A perturbation boundary element method (BEM) for steady-state flow problems in heterogeneous media is developed, and various aspects of the model are tested. The perturbation series obtained by the BEM gives the complete solution for the original governing equation. For highly heterogeneous problems. Pade approximants are effective to accelerate the rate of convergence and to convert a divergent series into a convergent solution. Because the model does not suffer from numerical dispersion and grid orientation effects, it is possible to track streamlines accurately.