In this study, we present an efficient approach, called the probabilistic collocation method (PCM), for uncertainty analysis of flow in unsaturated zones, in which the constitutive relationship between the pressure head and the unsaturated conductivity is assumed to follow the van Genuchten-Mualem model. Spatial variability of soil parameters leads to uncertainty in predicting flow behaviors. The aim is to quantify the uncertainty associated with flow quantities such as the pressure head and the effective saturation. In the proposed approach, input random fields, i.e., the soil parameters, are represented via the Karhunen-Loeve expansion, and the flow quantities are expressed by polynomial chaos expansions (PCEs). The coefficients in the PCEs are determined by solving the equations for a set of carefully selected collocation points in the probability space. To illustrate this approach, we use two-dimensional examples with different input variances and correlation scales and under steady state and transient conditions. We also demonstrate how to deal with multiple-input random parameters. To validate the PCM, we compare the resulting mean and variance of the flow quantities with those from Monte Carlo (MC) simulations. The comparison reveals that the PCM can accurately estimate the flow statistics with a much smaller computational effort than the MC.