These models have three mobile phases and at least three components. In three dimensions, Huyakorn et al. [1994] developed a three-phase three-component dual-porosity model suitable for fractured or non-fractured media. There was one NAPL component and no mass transfer, which avoided the burden of compositional calculations but limited the applicability of the model. The primary variables were one pressure and two saturations, as in typical black-oil petroleum reservoir models. A fully implicit mass-lumped upstream-weighted finite-element discretization with hexahedral bricks allowed for 7-, 11-, and 27-point connectivities and implicit wells. The implicit equations were solved by Newton-Raphson iteration, treating linear algebraic equations with a block-Orthomin procedure preconditioned by incomplete factorization (ILU). Positive transmissivities were assured by the formulation. Building on an earlier model of Forsyth and Shao [1991], Forsyth [1993] presented a two-dimensional three-phase four-component model to handle two contaminant components and steam. NAPL could flow in all three phases, water in the water and gas phases, and air in the gas phase only, with equilibrium mass transfer governed by empirical fluid-property correlations. The gas phase was assumed to be always present, with a small fictitious saturation if necessary, and on this basis variable substitution was implemented to treat phase appearance and disappearance. This forces the substitution logic into the flow and transport equations, unlike generalized compositional petroleum reservoir models that confine this logic to a fluid-property module by using total component concentrations as primary variables. A fully implicit (except for some dispersion terms) mass-lumped upstream-weighted control-volume finite-element method generated nonlinear equations that were solved by Newton-Raphson iteration with an ILU-preconditioned conjugate-gradient algorithm for linear equations. A major point of the paper was a procedure for modifying some dispersion coefficients in order to avoid negative transmissivities and obtain a maximum principle, which prevents nonphysical oscillations.