Inverse methods are increasingly used to estimate the hydraulic properties of unsaturated soils. The method generally uses a weighted least-squares approach in which numerically simulated data are fitted to measured data. In this study we used inverse methods to estimate the unsaturated soil hydraulic properties from continuous and multistep column outflow experiments. The method employs piecewise polynomial functions to obtain a free-form parameterization of the hydraulic properties, rather than fixed functional forms typical of the van Genuchten–Mualem and Brooks–Corey–Burdine models. For the polynomial functions we used quadratic B-splines and piecewise cubic Hermite interpolation. The method leads to local parameterizations that can also be hierarchic, depending on the invoked number of degrees of freedom. Since a suitable number of degrees of freedom cannot be defined a priori, we embedded the estimation method into a multilevel procedure, which also included a stability analysis, based on singular value decomposition of the sensitivity matrix. The optimization procedure was made more stable by imposing monotonicity constraints on the hydraulic functions. Tests with synthetic and measured data from column outflow experiments show the validity and robustness of the method.