Robust and detailed quantitative prediction of soil organic carbon (SOC) is of great significance to studying the carbon budget, soil management and decision-making. Spatial variations of SOC content were modelled using 863 soil profiles and a set of 22 environmental covariates representing relief, bioclimate variables and remote sensing data. The article provided the results of 3D modeling of SOC content in several soil layers (0–5, 5–15, 15–30, 30–60 and 60–100 cm) for the territory of the Russian Federation with 500 m spatial resolution. Machine learning framework was used, with random forest and spatial cross-validation techniques (150 km blocks) to handle the spatial autocorrelation of the training points. Compared with randomized cross-validation (R2 0.66, Concordance Correlation Coefficient (CCC) 0.79, RMSE 0.99 g/kg), using spatial cross-validation to predict the SOC content yielded less accurate results — R2 0.45, CCC 0.63, RMSE 1.41 g/kg. Regarding the importance of the variables, soil depth and temperature seasonality were major contributors to the SOC content prediction, followed by the EVI, 7 (MIR) MODIS band, and the topographic wetness index. The model was next evaluated with procedure so-called „area of applicability“ (AOA) of prediction model — the areas for which we cannot estimate prediction quality. AOA spatial distribution showed that the feature space not represented by training data is located in the mountain provinces. The proposed framework can be used for SOC modeling with a limited soil profile number, and it is provides a reproducible approach for long-term SOC monitoring.