Only incomplete knowledge of soil hydraulic properties is usually available since detailed sampling is rarely carried out in most soil surveys due to various constraints. Geostatistical interpolation techniques can then be used to estimate soil hydraulic properties at unsampled locations from existing observations. In addition to saturated hydraulic conductivities needed for water flow in saturated systems, knowledge of functional relationships describing soil hydraulic properties, i.e., retention and conductivity functions, are needed for unsaturated water flow. The soil water retention curve is usually parameterized using a simple functional form. These analytical models are convenient as they significantly reduce the number of variables (or parameters), compared to the number of measurements used to construct the original retention curve. The performance of two geostatistical approaches, parametric (P) and non-parametric (NP), are compared in this study to estimate parameters of the water retention function at unsampled locations. The impact of predicted parameters on simulations of variably-saturated water flow is also investigated. Water retention data used in this study were reported by Wierenga et al. (1989) for the Las Cruces trench site experiment, in which 450 soil samples were taken from nine layers (50 equally spaced samples per layer) from a 24.6-m by 6.0-m trench wall. Results show that for most water retention models the NP approach outperforms the P approach in terms of prediction errors. Unsaturated water flow was then simulated using HYDRUS-2D on kriged fields. Although both approaches show similar infiltration patterns, outliers can greatly impact predicted flow paths. Results show that the commonly and conveniently used parametric approach to geostatistically estimate water retention curves has to be used with care, while the NP approach leads systematically to acceptable results.