We propose a new method for the solution of stochastic unsaturated flow problems in randomly heterogeneous soils which avoids linearizing the governing flow equations or the soil constitutive relations, and places no theoretical limit on the variance of constitutive parameters. The proposed method applies to a broad class of soils with flow properties that scale according to a linearly separable model provided the dimensionless pressure head psi has a near-Gaussian distribution. Upon treating psi as a multivariate Gaussian function, we obtain a closed system of coupled nonlinear differential equations for the first and second moments of pressure head. We apply this Gaussian closure to steady-state unsaturated flow through a randomly stratified soil with hydraulic conductivity that varies exponentially with alpha psi where psi =(1/alpha)psi is dimensional pressure head and alpha is a random field with given statistical properties. In one-dimensional media, we obtain good agreement between Gaussian closure and Monte Carlo results for the mean and variance of psi over a wide range of parameters provided that the spatial variability of alpha is small. We then provide an outline of how the technique can be extended to two- and three-dimensional flow domains. Our solution provides considerable insight into the analytical behavior of the stochastic flow problem.