Unsaturated steady flow in three dimensional bounded heterogeneous formations is modelled by considering the hydraulic parameters of the conductivity curve as stationary Random Space Functions (RSFs). The corresponding stochastic governing equations are solved by means of a perturbation approach which regards the variances of the logarithms of the saturated hydraulic conductivity (Ks), and the pore scale distribution of the Gardner (1958) model as small quantities. Analytical closed forms (requiring only few quadratures) for the spatial moments of the pressure head are derived, without invoking the unit mean hydraulic assumption. It is found that the flow variables are nonstationary near the boundary (i.e. the water table), and approach stationarity as the vertical distance from the water table increases. The stationary limits and the critical vertical distance at which stationarity is attained depend on soil types as well as infiltrating rates. The smaller the rate is, the larger the critical distance; and the coarser the soil texture is, the smaller the distance. Our results generalize the study of Indelman et al. (1993), and provide a general methodology to tackle more complex unsaturated flows