A three-dimension boundary element method is developed to study unsaturated groundwater flows in heterogeneous porous media. The method is based on the Kirchhoff transformation of variables which transforms the nonlinear Richards equation into the Fokker-Planck equation. Several benchmark problems are considered involving cavity infiltration and exclusion to demonstrate the reliability of the boundary element formulation. As an example application problem, the effective conductivity of a porous medium containing randomly positioned impermeable spherical inclusions is considered. ¿ American Geophysical Union 1994