A physically motivated model is indispensable for a successful analysis of the impact of leaching from nuclear waste storage sites on the environment and public health. While most analyses use the single-phase flow model for modelling unsaturated flow and solute transport, the two-phase flow model considering the resistance of gas to water flow is a more realistic one. The effect of the two-phase flow model on the water content is theoretically investigated first in this study. Then, by combining a geostatistical generator using the turning bands method and a multi-phase transport code TOUGH2, an automatic process is used for Monte Carlo simulation of the solute transport. This stochastic approach is applied to a potentially polluted site by low-level nuclear waste in Taiwan. In the simulation, the saturated hydraulic conductivity is treated as the random variable. The stochastic approach provides a probabilistic assessment of contamination. The results show that even though water content from the two-phase flow model is only 1.5% less than the one from the single-phase flow model, the two-phase flow causes a slower movement but a wider lateral spreading of the plume in the unsaturated zone. The stochastic approach provides useful probability information which is not available from the deterministic approach. The probability assessment of groundwater contamination provides the basis for more informed waste management, better environmental assessment and improved evaluation of impact on public health.