Phreatic flow in heterogeneous aquifers is analyzed by treating hydraulic conductivity as a random field with known statistics. A set of equations for the first and second ensemble moments of hydraulic head and phreatic surface is derived. These equations allow one to predict the behavior of phreatic aquifers, as well as to assess the uncertainty associated with such predictions. Perturbation analysis in variance of log hydraulic conductivity is employed to close the moments equations. This leads to a recursive initial boundary value problem.