Soil hydraulic properties (SHPs) are commonly determined in soil samples with replicas. Whether these replicas are taken at a same location to represent a specific point or at several locations to represent a larger area, the results should be merged into a final data set to be used in modeling. For this data set to be representative, standard errors and a correlation matrix must be considered in the merging process. We present a method to perform this merging and give an example using stochastic realizations of van Genuchten-Mualem (VGM) parameters generated by Cholesky decomposition to merge the SHP and associated statistics into a merged parameter set. To do so, we used VGM parameters obtained at sample scale in three replicas from a Brazilian savanna soil through inverse modeling of laboratory evaporation experiments. The effectiveness and representativeness of the proposed methodology were evaluated by observing the frequency distribution of different levels of output, comparing individual and merged sample properties. The outputs include VGM parameters, retention and conductivity characteristics, and water balance components stochastically predicted by a hydrological model. The performed stochastic merging correctly represented the variability of the combined replicas, especially with respect to hydrological model outputs of soil water balance components. Using the mean hydraulic property parameter values to deterministically predict water balance components may yield values that are substantially different from the mean values of stochastic realizations. This suggests that the deterministic prediction using mean parameter values in vadose zone hydrological modeling may result in unrepresentative outputs.