Numerical modeling of the migration of three-phase immiscible fluid flow in variably saturated zones is challenging due to the different behavior of the system between unsaturated and saturated zones. This behavior results in the use of different numerical methods for the numerical simulation of the fluid flow depending on whether it is in the unsaturated or saturated zones. This paper shows that using a high-resolution shock-capturing conservative method to resolve the nonlinear governing coupled partial differential equations of a three-phase immiscible fluid flow allows the numerical simulation of the system through both zones providing a unitary vision (and resolution) of the migration of an immiscible contaminant problem within a porous medium. In particular, using different initial scenarios (including impermeable "lenses" in heterogeneous aquifers), three-dimensional numerical simulation results are presented on the temporal evolution of the contaminant migration following the saturation profiles of the three-phases fluids flow in variably saturated zones. It is considered either light nonaqueous phase liquid with a density less than the water, or dense nonaqueous phase liquid, which has densities greater than the water initially released in unsaturated dry soil. Our study shows that the fate of the migration of immiscible contaminants in variably saturated zones can be accurately described, using a unique mathematical conservative model, with different evolution depending on the value of the system's physical parameters, including the contaminant density, and accurately tracking the evolution of the sharp (shock) contaminant front.