The present work describes a model developed for multiphase transport in the subsurface. The system under consideration comprises three phases including the immobile solid phase composed of soil grains, the aqueous phase flowing through the bed of soil grains, and the non-aqueous phase of liquids (NAPLs) that may consist of several mutually soluble compounds. Numerous field data and experimental results have indicated that the residual nonwetting fluids of the NAPLs in groundwater systems are trapped, i.e., completely surrounded by the wetting aqueous phase. In the model, therefore, the NAPLs are treated as discrete blobs with generalized local-size distribution, while the aqueous phase is assumed to be a continuum. In addition, the model takes into account the surface area-to-volume ratio and the ratio of the aqueous contacting area to the overall surface area of NAPL blobs. Interactions between the two liquid phases manifest themselves in the governing equations of the aqueous phase as an integral over all NAPL blobs. The rate-controlling factors in the transport process have been analyzed. By resorting to a semi-implicit finite-difference algorithm, numerical studies have been carried out to examine characteristics of the system, such as time-dependent concentration profiles, remaining quantities of the contaminant, and variations of the NAPL blob size. The results of simulation reveal that the initial blob size of the NAPLs and hydrodynamic conditions profoundly affect the rates of dissolution and transport; desorption and dissolution occur simultaneously, but the dissolution is completed first; and the NAPL attenuates via a moving front in the soil bed under convection-dispersion control. The model predictions for toluene dissolution in a glass bead column agree well with the available experimental data.