The knowledge about the spatial distribution of soil organic carbon stock (SOCS) helps in sustainable land-use management and ecosystem functioning. No such study has been attempted in the complex topography and land use of Himalayas, which is associated with great spatial heterogeneity and uncertainties. Therefore, in this study digital soil mapping (DSM) was used to predict and evaluate the spatial distribution of SOCS using advanced geostatistical methods and a machine learning algorithm in the Himalayan region of Jammu and Kashmir, India. Eighty-three soil samples were collected across different land uses. Auxiliary variables (spectral indices and topographic parameters) derived from satellite data were used as predictors. Geostatistical methods—ordinary kriging (OK) and regression kriging (RK)—and a machine learning method—random forest (RF)—were used for assessing the spatial distribution and variability of SOCS with inter-comparison of models for their prediction performance. The best fit model validation criteria used were coefficient of determination (R2) and root mean square error (RMSE) with resulting maps validated by cross-validation. The SOCS concentration varied from 1.12 Mg/ha to 70.60 Mg/ha. The semivariogram analysis of OK and RK indicated moderate spatial dependence. RF (RMSE = 8.21) performed better than OK (RMSE = 15.60) and RK (RMSE = 17.73) while OK performed better than RK. Therefore, it may be concluded that RF provides better estimation and spatial variability of SOCS; however, further selection and choice of auxiliary variables and higher soil sampling density could improve the accuracy of RK prediction.