Water flow analyses under transient soil hydraulic conditions require knowledge of the soil hydraulic properties. These constitutive relationships, named soil-water characteristic curve (SWCC) and hydraulic conductivity function (HCF) are described through empirical models which generally have several parameters that must be calibrated against collected data. Many of the parameters in SWCC and HCF models cannot be directly measured in field or laboratory but can only be meaningfully inferred from collected data and inverse modeling. In order to obtain the soil parameters with the inverse process, a local or global optimization algorithm may be applied. Global optimizations are more capable of fiding optimum parameters, however the direct solution through numerical modeling are time consuming. Therefore, analytical solutions (surrogate models) may overcome this shortcomming by accelerating the optimization process. In this work we introduce Evolutionary Polynomial Regression (EPR) as a tool to develop surrogate models of the physically-based unsaturated flow. A rich dataset of soil hydraulic parameters is used to calibrate our surrogate model, and real-world data are then utilized to validate our methodology. Our results demonstrate that the EPR model predicts accurately the observed pressure head data. The model simulations are shown to be in good agreement with the Hydrus software package.