Many spatial interpolation methods perform well for gentle terrains when producing spatially continuous surfaces based on ground point data. However, few interpolation methods perform satisfactorily for complex terrains. Our objective in the present study was to analyze the suitability of several popular interpolation methods for complex terrains and propose an optimal method. A data set of 153 soil water profiles (1 m) from the semiarid hilly gully Loess Plateau of China was used, generated under a wide range of land use types, vegetation types and topographic positions. Four spatial interpolation methods, including ordinary kriging, inverse distance weighting, linear regression and regression kriging were used for modeling, randomly partitioning the data set into 2/3 for model fit and 1/3 for independent testing. The performance of each method was assessed quantitatively in terms of mean-absolute-percentage-error, root-mean-square-error, and goodness-of-prediction statistic. The results showed that the prediction accuracy differed significantly between each method in complex terrain. The ordinary kriging and inverse distance weighted methods performed poorly due to the poor spatial autocorrelation of soil moisture at small catchment scale with complex terrain, where the environmental impact factors were discontinuous in space. The linear regression model was much more suitable to the complex terrain than the former two distance-based methods, but the predicted soil moisture changed too sharply near the boundary of the land use types and junction of the sunny (southern) and shady (northern) slopes, which was inconsistent with reality because soil moisture should change gradually in short distance due to its mobility in soil. The most optimal interpolation method in this study for the complex terrain was the hybrid regression kriging, which produced a detailed, reasonable prediction map with better accuracy and prediction effectiveness.