Detailed maps of regional spatial distribution of soil organic carbon (SOC) are needed to guide sustainable soil uses and management decisions. Interpolation methods based on spatial auto-correlations, environmental covariates, or hybrid methods are commonly used to predict SOC maps. Many of these methods perform well for gentle terrains. However, it is unknown how these methods perform to capture SOC variations in complex terrains, especially areas of which land uses are interrupted by human activities, such as the Loess Plateau of China. This study compared four interpolations or predictive methods including ordinary kriging (OK), regression kriging, ordinary kriging integrated with land-use type (OK_LU) and a soil land inference model (SoLIM). The purpose of this study is to find appropriate methods, which are suitable to the complex terrain in Loess Plateau region of China. The study area was a typical watershed in Loess Plateau with complex hilly–gully terrain and various land-use types. A field sampling dataset of 200 points was partitioned into 1/2 for model building and 1/2 for accuracy validation in a random way. Nine environmental covariates were selected: land-use types, digital elevation model, solar radiation, slope degree, slope aspect, plan curvature, profile curvature, surface area ratio, and topographic wetness index. The mean absolute percentage error, root mean square error, and goodness-of-prediction statistic value were selected to evaluate mapping results. The results showed that the use of easily obtained environmental covariates, land-use types and terrain variables improved accuracies of SOC interpolation, which will be of interests for related research of similar environments in the Loess Plateau. SoLIM and OK_LU can be two suitable and efficient methods, which produced detailed, reasonable maps with higher accuracy and prediction effectiveness, for the study area and similar areas in the Loess Plateau.