At sites where a dense nonaqueous phase liquid (DNAPL) was spilled or released into the subsurface, estimates of the mass of DNAPL contained in the subsurface from core or monitoring well data, either in the nonaqueous or aqueous phase, can be highly uncertain because of the erratic distribution of the DNAPL due to geologic heterogeneity. In this paper, a multiphase compositional model is applied to simulate, in detail, the DNAPL saturations and aqueous-phase plume migration in a highly characterized, heterogeneous glaciofluvial aquifer, the permeability and porosity data of which were collected by researchers at the University of Tübingen, Germany. The DNAPL saturation distribution and the aqueous-phase contaminant mole fractions are then reconstructed by sampling the data from the forward simulation results using two alternate approaches, each with different degrees of sampling conditioning. To reconstruct the DNAPL source zone architecture, the aqueous-phase plume configuration, and the contaminant mass in each phase, one method employs the novel transition probability/Markov chain approach (TP/MC), while the other involves a traditional variogram analysis of the sampled data followed by ordinary kriging. The TP/MC method is typically used for facies and/or hydraulic conductivity reconstruction, but here we explore the applicability of the TP/MC method for the reconstruction of DNAPL source zones and aqueous-phase plumes. The reconstructed geometry of the DNAPL source zone, the dissolved contaminant plume, and the estimated mass in each phase are compared using the two different geostatistical modeling approaches and for various degrees of data sampling from the results of the forward simulation. It is demonstrated that the TP/MC modeling technique is robust and accurate and is a preferable alternative compared to ordinary kriging for the reconstruction of DNAPL saturation patterns and dissolved-phase contaminant plumes.