The Karhunen-Loeve moment equation (KLME) approach is implemented to model stochastic transient water-NAPL two-phase flow in heterogeneous subsurface media with random soil properties. To describe the constitutive relationships between water saturation, capillary pressure, and phase relative permeability, the widely used van Genuchten model and Parker and Lenhard models are adopted. The log-transformed intrinsic permeability, soil pore size distribution, and van Genuchten fitting parameter n are treated as normally distributed stochastic variables with a separable exponential covariance model. The perturbation part of these three log-transformed variables is decomposed via Karhunen-Loeve expansion. The dependent variables (phase pressure, phase mobility, and capillary pressure) are expanded by polynomial expansions and the perturbation method. Incorporating these expansions of random soil properties variables and dependent variables into the governing equations yields a series of differential equations in different orders. We construct the moments of the dependent variables from the solutions of these differential equations. We demonstrate the stochastic model with two-dimensional examples of transient two-phase flow. We also conduct Monte Carlo simulations using the finite element heat and mass (FEHM) transfer code, whose results are considered "true" solutions. The match between the results from FEHM and KLME indicates the validity of the proposed KLME application in transient two-phase flow. The computational efficiency of the KLME approach over Monte Carlo methods is at least an order of magnitude for transient two-phase flow problems.