Simulation of three-dimensional flow of immiscible fluids within and below the unsaturated zone

This paper presents a two‐phase flow model based on a three‐dimensional, finite‐difference formulation. As three‐dimensional simulations can require substantial computer effort, a numerical technique that takes advantage of vector and parallel processing computer architecture is developed. The model is posed in terms of water saturation and nonwetting fluid pressure. It uses three‐phase capillary pressure and relative permeability relationships to permit simulation within or below the unsaturated zone. A modified formulation of slice successive overtaxation (an iterative matrix solution technique) is introduced. This technique is designed to use parallel processing capabilities of new computers. The model is applied to immiscible fluid flow at two chemical waste landfills near Niagara Falls, New York. At both sites, denser than water, nonaqueous liquids (NAPLs) are present in the groundwater regimes in relatively large quantities. The model applications address several technical concerns at the two sites, including the effectiveness of clay as a geologic barrier to NAPL migration owing to capillary pressure forces, the three‐dimensional aspects of dense NAPL flow, and the sensitivity of NAPL recovery in pumping wells due to various hydrogeologic and fluid properties. The results of the applications show that (1) even under a downward hydraulic gradient, natural differences in capillary pressure relationships for different lithologies can prevent downward migration of NAPL, (2) without any lithologic‐capillary barrier, an upward hydraulic gradient induced by a de watering system can prevent downward migration of NAPL, (3) NAPL recovery at wells is sensitive to relative permeability, a relationship that requires field calibration in many settings, and (4) the three‐dimensional aspects of two‐phase flow and hydrogeologic heterogeneity require explicit treatment in many settings.