Saturated hydraulic conductivity (Ks) is one of the key parameters in modeling solute and water movement in the vadose zone. Field and laboratory measurement of Ks is time consuming, and hence is not practical for characterizing the large spatial and temporal variability of Ks As an alternative to direct measurements, pedotransfer functions (PTFs), which estimate Ks from readily available soil data, are being widely adopted. This study explores the utility of a promising data‐driven method, namely, genetic programming (GP), to develop PTFs for estimating Ks from sand, silt, and clay contents and bulk density (Db). A data set from the Unsaturated Soil Hydraulic Database (UNSODA) was considered in this study. The performance of the GP models were compared with the neural networks (NNs) model, as it is the most widely adopted method for developing PTFs. The uncertainty of the PTFs was evaluated by combining the GP and the NN models, using the nonparametric bootstrap method. Results from the study indicate that GP appears to be a promising tool for developing PTFs for estimating Ks The better performance of the GP model may be attributed to the ability of GP to optimize both the model structure and its parameters in unison. For the PTFs developed using GP, the uncertainty due to model structure is shown to be more than the uncertainty due to model parameters. Moreover, the results indicate that it is difficult, if not impossible, to achieve better prediction and less uncertainty simultaneously.