Soil heterogeneity and the lack of detailed site characterization are two ubiquitous factors that render predictions of flow and transport in the vadose zone inherently uncertain. We employ the Green‐Ampt model of infiltration and the Dagan‐ Bresler statistical parameterization of soil properties to compute probability density functions (PDFs) of infiltration rate and infiltration depth. By going beyond uncertainty quantification approaches based on mean and variance of system states, these PDF solutions enable one to evaluate probabilities of rare events that are required for probabilistic risk assessment. We investigate the temporal evolution of the PDFs of infiltration depth and corresponding infiltration rate, the relative importance of uncertainty in various hydraulic parameters and their cross-correlation, and the impact of the choice of a functional form of the hydraulic function. Soil heterogeneity and the lack of detailed site characterization are two ubiquitous factors that hamper one’s ability to predict flow and transport in the vadose zone. The continuing progress in data acquisition notwithstanding, measurements of hydraulic properties of partially saturated media remain scarce and prone to measurement and interpretive errors. Consequently, spatial distributions of hydraulic parameters (saturated and relative hydraulic conductivities, and parameters in retention curves) are typically uncertain and their statistical properties are subject to considerable debate. Despite some reservations, e.g., [1, 2], it has become common to treat saturated hydraulic conductivityK s (x)as a multivariate log-normal random field whose ensemble statistics (e.g., mean, variance, and correlation length) can be inferred from spatially distributed data by means of geostatistics. No such consensus exists about statistical distributions of various hydraulic parameters entering relative hydraulic conductivity and retention curves. For example, various data analyses concluded that spatial variability of a soil parameterα G (x)in the Gardner model of relative conductivity, which is often referred to as the reciprocal of the macroscopic capillary length, exhibits either a normal [3] or log-normal [4] distribution and is either correlated [5] or uncorrelated [3] withK s . We defer a more detailed review of the statistical properties of bothα G (x)and parameters in the van Genuchten model of relative conductivity until Section 2. Here, it suffices to say that any approach to uncertainty quantification for flow and transport in the vadose zone must be flexible enough to accommodate arbitrary statistical distributions of soil properties. Statistical treatment of hydraulic parameters renders corresponding flow equation stochastic. Solutions of these equations are probability density functions (PDFs) of system states (water content, pressure, and macroscopic flow velocity) and can be used not only to predictflowin heterogeneous partially saturated porous media but also to quantify