We consider free surface flow in random porous media by treating hydraulic conductivity of a medium as a random field with known statistics. We start by recasting the boundary-value problem in the form of an integral equation where the parameters and domain of integration are random. Our analysis of this equation consists of expanding the random integrals in Taylor's series about the mean position of the free boundary and taking the ensemble mean. To quantify the uncertainty associated with such predictions, we also develop a set of integro-differential equations satisfied by the corresponding second ensemble moments. The resulting moment equations require closure approximations to be workable. We derive such closures by means of perturbation expansions in powers of the variance of the logarithm of hydraulic conductivity. Though this formally limits our solutions to mildly heterogeneous porous media, our analytical solutions for one-dimensional flows demonstrate that such perturbation expansions may remain robust for relatively large values of the variance of the logarithm of hydraulic conductivity.