We consider flow in partially saturated heterogeneous porous media with uncertain hydraulic parameters. By treating the saturated conductivity of the medium as a random field, we derive a set of deterministic equations for the statistics (ensemble mean and variance) of fluid pressure. This is done for three constitutive models that describe the nonlinear dependence of relative conductivity on pressure. We use the Kirchhoff transform to map Richards equation into a linear PDE and explore alternative closures for the resulting moment equations. Regardless of the type of nonlinearity, closure by perturbation is more accurate than closure based on the non-perturbative Gaussian mapping. We also demonstrate that predictability of unsaturated flow in heterogeneous porous media is enhanced by choosing either the Brooks–Corey or van Genuchten constitutive model over the Gardner model.