We consider steady state unsaturated flow in bounded randomly heterogeneous soils under the influence of random forcing terms. Our aim is to predict pressure heads and fluxes without resorting to Monte Carlo simulation, upscaling or linearization of the constitutive relationship between unsaturated hydraulic conductivity and pressure head. We represent this relationship through Gardner's exponential model, treating its exponent alpha as a random constant and saturated hydraulic conductivity, K_s, as a spatially correlated random field. This allows us to linearize the steady state unsaturated flow equations by means of the Kirchhoff transformation, integrate them in probability space, and obtain exact integro-differential equations for the conditional mean and variance-covariance of transformed pressure head and flux. We solve the latter for flow in the vertical plane, with a point source, by finite elements to second-order of approximation. Our solution compares favorably with conditional Monte Carlo simulations, even for soils that are strongly heterogeneous.