In evaluating risks associated with chemical spills on the ground surface or leaks from underground storage tanks, it is required to know the extent and degree of contamination in the subsurface. Under certain conditions the contamination front can be treated as sharp and quasi-steady during both initial migration away from a spill and during removal via in-situ extractive processes. The boundary element method (BEM) is a viable choice of numerical method of solution in such situations where only the movement of the boundary is required. BEM is used to solve the set of partial differential equations governing the subsurface flow of NAPL bound by a sharp surface or interface in the vertical plane. This is followed by moving the free boundary using Darcy's law. Unlike other methods the node movement is not restricted to the vertical direction. In order to obtain a high accuracy solution a node redistribution algorithm developed by Carey and Kennon  was implemented. It redistributed the nodes according to curvature of the free boundary while retaining the shape. Aquifers with homogeneous as well as piecewise homegeneous properties are considered. Results are presented for simulation of napl infiltration and spreading on an impermeable surface that agree well with laboratory infiltration experiments. Model parameters were obtained by separate measurements and were not fitted to the experimental data.