A numerical model has been developed to provide areal analyses of the three-dimensional spreading of immiscible liquids in groundwater systems. The model is intended for specific use in prediction of hydrocarbon spreading from subsurface leaks and spills and in design and performance evaluation of remedial schemes. The mathematical formulation is based on vertical integration of the three-dimensional two-phase flow equations and incorporation of the concept of gravity-capillary vertical equilibrium (GCVE) in which a vertical balance of gravitational and capillary forces is assumed. History-dependent pseudo functions of capillary pressure and relative permeabilities are introduced for the present GCVE model. A direct numerical procedure for evaluating the pseudo functions is derived and verified using an example given in the petroleum literature. An overview of numerical techniques for solving the nonlinear governing equations is presented. Special schemes for handling production wells are derived. Simulation examples are provided to verify, validate, and demonstrate utility of the model. Numerical results obtained from the GCVE model are compared with analytical and rigorous multiphase three-dimensional numerical solutions. The model is further validated using data from a laboratory investigation on waterflood, five-spot well performance. The results of these comparisons show the validity of the GCVE modeling assumptions and accuracy and robustness of the proposed formulations and computational schemes. The numerical study also indicates that the present model is highly efficient and suitable for field simulations on personal computers or workstations.