A stochastic model for saturated-unsaturated flow is developed based on the combination of the Karhunen-Loeve expansion of the input random soil properties with a perturbation method. The saturated hydraulic conductivity ks(x) is assumed to be log-normal random functions, expressed by f(x). f(x) is decomposed as infinite series in a set of orthogonal normal random variables by the Karhunen-Loeve (KL) expansion and the pressure head is expand as polynomial chaos with the same set of orthogonal random variables. With these expansions, the stochastic saturated-unsaturated flow equation and the corresponding initial and boundary conditions are represented by a series of deterministic partial differential equations which can be solved subsequently by a suitable numerical method. Some examples are given to show the reliability and efficiency of the proposed method.