Spatial heterogeneity of geologic media leads to uncertainty in predicting both flow and transport in the vadose zone. In this work an efficient and flexible, combined analyticalnumerical Monte Carlo approach is developed for the analysis of steady-state flow and transient transport processes in highly heterogeneous, variably saturated porous media. The approach is also used for the investigation of the validity of linear, first order analytical stochastic models. With the Monte Carlo analysis accurate estimates of the ensemble conductivity, head, velocity, and concentration mean and covariance are obtained; the statistical moments describing displacement of solute plumes, solute breakthrough at a compliance surface, and time of first exceedance of a given solute flux level are analyzed; and the cumulative probability density functions for solute flux across a compliance surface are investigated. The results of the Monte Carlo analysis show that for very heterogeneous flow fields, and particularly in anisotropie soils, The linearized, analytical predictions of soil water tension and soil moisture flux become erroneous. Analytical, linearized Lagrangian transport models also overestimate both the longitudinal and the transverse spreading of the mean solute plume in very heterogeneous soils and in dry soils. A combined analytical-numerical conditional simulation algorithm is developed to estimate the impact of in-situ soil hydraulic measurements on reducing the uncertainty of concentration and solute flux predictions. In soils with large spatial variability and in dry soils, soil water tension measurements significantly reduce the uncertainty in the predicted solute concentration. Saturated hydraulic conductivity data are valuable in relatively wet soils. A combination of tension and saturated hydraulic conductivity data gives the best results, especially if some data are available on the unsaturated hydraulic conductivity function. It is also found that if soil heterogeneity is large, the conditional spatial moments of inertia of the mean concentration plume and the conditional mean breakthrough curves are poor means of depicting the actual solute plume distribution and the actual solute flux. Nevertheless, conditional simulation is one of the most rational approaches for modeling unsaturated flow and transport, if in-situ data are available.