A stochastic unsaturated water flow model is developed for heterogeneous soils subject to a transient flow regime. Equations are developed for a fully three-dimensional soil profile, and results are presented for an example one-dimensional problem. The model predicts the mean and covariance of the soil water content, Darcy's flux, and pore water velocity as a function of the boundary flux and saturated hydraulic conductivity statistics. The statistics of the pore water velocity can be used to predict solute transport in soils, as shown by Foussereau et al. . Approximate flow-related moment equations are solved analytically in the Laplace domain. Then, the analytical results are numerically inverted using a modified fast Fourier transform algorithm. The model predictions are compared to results obtained from Monte Carlo simulations for two different boundary flux patterns characteristic of humid climates and two different soil types (a fine sand and a sandy loam). Comparing the approximate solutions of the statistical moments to the outputs of the Monte Carlo simulations shows (1) the dominance of the boundary flux variability over that of the saturated conductivity on the overall prediction uncertainty, particularly at shallow depths, and (2) the good performance of the stochastic unsaturated flow model, particularly for fine-textured soils subject to boundary fluxes with coefficients of variation up to ~1.5. As the boundary flux coefficient of variation increases and the soil becomes coarser, the model performance deteriorates because the flow system becomes significantly more nonlinear. ¿ 2000 American Geophysical Union