【Objective】Soil organic carbon (SOC) is an important indicator of soil fertility and plays a fundamental role in the terrestrial ecosystem carbon cycle. As one of the primary environmental factors in digital soil mapping (DSM), landform elements are irreplaceable in predicting SOC. The purpose of this study was to simulate the complex and nonlinear relationship between SOC and environmental variables and evaluate the importance of each variable to accuracy in SOC mapping. 【Method】We applied machine learning techniques to map SOC content in a small watershed (1: 25000) of Huangmei Town, Jurong City using high-resolution landform elements classification maps known as geomorphons, digital elevation model (DEM) derivatives, optical and synthetic aperture radar (SAR) remote sensing data. The performance of all geomorphon (GM) variables under different hyperparameter settings was evaluated to predict SOC content. Three machine-learners including bagged classification and regression tree (Bagged CART), random forest (RF) and Cubist were used to construct predictive models of SOC content based on 74 soil samples and different combinations of environmental covariates. Model A, Model B, and Model D included only GM variables, DEM derivatives, and remote sensing variables, respectively. Model B was a combination of GM data and DEM derivatives, while Model E included all predictor variables. The performance of these models was evaluated based on a 10-fold cross-validation method by four statistical indicators. Concordance index (C‑index), root mean square errors (RMSE), bias and coefficient of determination (R²) of the three models were worked out for evaluation of the accuracy of their predictions. The best model was screening-out for mapping SOC in the study area based on the raster datasets of all environmental variables. 【Result】Overall, the Cubist model performed better than RF and Bagged CART, and these models yielded similar spatial distribution patterns of SOC, i.e. an ascending trend from the northern hilly area to the southern flatter land of the study area. Our results showed that more accurate predictions of SOC content were provided with the introduction of GM variables than individual DEM derivatives. The GM map with 20 cells search radius (L) and 5° flatness threshold (t) showed the highest relative importance within four GM variables in three models. The Cubist‑E model that functioned based on GM landform elements classification variables, DEM derivatives and remote sensing variables was much better than the others in performance and could explain most of the spatial heterogeneity of SOC (R²= 0.53). Also, the prediction accuracy changed with and without the GM predictors with the R² for estimating SOC content using the Cubist model increasing by 14.3%. The SOC contents of the hilly region predicted with the Cubist‑E model ranged from 5.65 to 13.31 g·kg^–1. In addition, topographic variables were the main explanatory variables for SOC predictions and the multi-resolution index of valley bottom flatness (MRVBF) and elevation were assigned as the two most important variables. 【Conclusion】The Cubist model that functions based on GM variables, DEM derivatives, as well as remote sensing variables, is a promising approach to predicting the spatial distribution of SOC in hilly regions at a small watershed scale. The results of this study illustrate the potential of GM landform elements classification data as input when developing SOC prediction models.