A numerical model for the analysis of uncertainty propagation in flow through unsaturated soils is developed. This model is based on a first-order Taylor series expansion of the discretized Richards equation. Soil hydrologic properties (the saturated hydraulic conductivity and the pore size distribution parameter) are assumed to be stochastic processes in space. The surface boundary conditions can be considered as deterministic variables in time or stochastic time series. Spectral analysis and Monte Carlo simulations were used to verify this numerical model for flow under both steady and transient conditions. The model is then used to examine the effect of uncertainty in boundary conditions and the effect of heterogeneity on the pressure head and flux variance profiles at various times for one-dimensional vertical flow cases. Dependence of pressure head variance on the flow conditions (drying or wetting) is examined. On the basis of the analysis it is found that the propagation of the head variance is similar to that of the concentration variance for solute transport in saturated aquifers. The head variance is proportional to the mean pressure gradient, and thus large head variances are associated with the wetting and the drying front of a moisture pulse. The peak head variance is smaller at the wetting front than it is at the drying front. This difference is attributed to the difference in the magnitude of mean hydraulic gradient and should not necessarily be interpreted as a hysteresis effect. In addition, it is shown how the variance of the flux of a moisture pulse increases with depth. ¿ 1999 American Geophysical Union