A newly developed multiphase mixture model is applied to numerically investigate infiltration and transport of nonaqueous phase liquids (NAPLs) in the unsaturated subsurface. Simultaneous flows of liquid and gas phases, solutal convection and the associated organic vapor transport in the gas phase are accounted for in the numerical model. It is shown that the numerical modeling of complex and interactive transport processes in multiphase, multicomponent systems can become computationally less intensive if based on this new model. The numerical results for three common contaminants reveal that the two-phase zone resulting from a NAPL spill is generally characterized by two distinct regions: an elongated pancake-shaped lens of high NAPL saturation which floats over the water table, and the area above it at lower saturation which was previously swept by the NAPL infiltration front. These predicted features are consistent with previous experimental observations for lighter than water NAPLs. It is also found that the vapor phase transport responsible for large-scale contamination is mainly caused by the displacement flow occurring during NAPL infiltration. In comparison, the density-driven flow due to evaporation of heavier organic chemicals is less important for systems with mobile NAPLs. Finally, several areas where future research is needed are discussed.