Water retention characteristic curve (WRCC) is a primary input for modeling the unsaturated behavior of geomaterials. The inadequacy of measurement methodologies in capturing the measurable range of the suction of geomaterials along with material variability induce uncertainties in the quantification of WRCC. The goal of this study is to demonstrate the usefulness of copulas for quantifying the uncertainties in WRCC attributed to the aforementioned reasons. Four instruments measuring different ranges of suction of four fly ashes were used to generate 54 combinations of measured data that would represent possible cases of data-induced uncertainties. van Genuchten WRCC equation parameters (a and n) were determined for these combinations. Statistical analysis of the parameters determined from the 54 combinations indicates a negative correlation with each other; therefore, they represent a bivariate random vector. Four commonly used copulas, Gaussian, Frank, Plackett, and No. 16, were utilized to construct the bivariate density function of the parameters, among which the Gaussian copula was found to be the best fit. The impact of the copula-based model developed in this study on two important unsaturated functions, the suction stress characteristic curve (SSCC) and hydraulic conductivity function (HCF) of fly ash, was investigated and found to be significant. The dispersion in SSCC ranges from monotonically increasing to sharp post peak reduction in the suction range of 20–50 kPa. The limiting state of suction required to reach the theoretical minimum permeability was found to vary widely in the range of 300–9,000 kPa. The study shows that a Gaussian copula with lognormal marginals can adequately capture the significant range of uncertainties associated with WRCC.