A formulation is presented for numerical modeling of light nonaqueous phase liquid (LNAPL) accumulation in observation wells and LNAPL removal from recovery wells. The formulation includes fluid gravity segregation, changing water and LNAPL levels in the well, and the corresponding changes in fluid storage in the well bore. The method was added to a three-dimensional finite difference model for hysteretic multiphase flow. The model predictions are compared to three-dimensional pilot scale experiments of LNAPL infiltration and recovery. Model predictions matched experimental results well, demonstrating the validity of the model formulation. Model results indicated that soil layering and lack of mechanical equilibrium had a significant influence on LNAPL thicknesses in wells. Thus empirical approaches based on assumptions of soil homogeneity and mechanical equilibrium would not have been appropriate for prediction of LNAPL volume from toluene thicknesses measured in wells. ¿ 2000 American Geophysical Union