Variability in the unsaturated soil hydraulic properties affects the behavior of water and solutes in the subsurface. When these properties are described by parametric functions, subsurface heterogeneity can be quantified by variations in the parameters of these functions. We propose a procedure to characterize systematically the statistical properties of each soil hydraulic parameter and the parameter correlations involved. First, the parameter with the best-defined probability distribution function (pdf) is identified. Curvilinear regression analysis identifies transformations that linearize the relationships between this reference parameter and the remaining parameters. From this we determine the joint pdf of the transformed parameters, which is then used to calculate the distributions of derived variates (functions of the hydraulic parameters). The procedure is applied to the hydraulic parameters of 140 samples in total taken from two layers of an agricultural soil profile. We developed analytical expressions for the pdfs of the derived variates (hydraulic conductivities and water contents at pressure heads between 0 and −5000 cm) from the multivariate parameter pdf. The numerical integration required to evaluate these expressions proved extremely cumbersome, thus reducing the robustness of the analytical expressions. We therefore performed a Monte Carlo simulation from which we determined the first two moments, as well as the skewness and excess of the derived variates. The estimated mean and standard deviation of the derived variates generally agreed well with values determined directly from the soil samples. Skewness and excess were estimated less accurately. Both the derived analytical pdfs and the Monte Carlo simulation showed markedly nonsymmetrical pdfs, suggesting that it is generally insufficient to limit statistical analyses to the first two moments only.