Understanding the vertical and lateral distribution of soil organic carbon (SOC) and soil organic carbon density (SOCD) is indispensable for soil use and environmental management because of their vital role in soil quality assessments. Primarily, they are needed in calculating soil organic carbon storage (SOCS). The objective of this research was to provide digital maps of SOC and SOCD variation as well as their uncertainties at multiple standardized depths (H1: 0–5, H2: 5–15, H3: 15–30, H4: 30–60 and H5: 60–100 cm) using a parsimonious model with optimized terrain‐related attributes and satellite‐derived data. SOCS were evaluated at soil subgroup levels. An area of about 808 km2 with varying elevation, plant cover and lithology from the Miandoab region, West Azerbaijan Province, Iran was selected as a case study area. A total of 386 soil samples were collected from 104 profiles comprising various soil genetic horizons. A continuous spline function was then fitted to the target properties in advance of creating a dataset at five standard depth intervals (following the GlobalSoilMap project). These were then grouped into three classes including top (H1), middle (H2, H3 and H4) and bottom (H5) depths to ease interpretation. Static and dynamic covariates (30‐m resolution) were derived from a digital elevation model (DEM) and a suite of Landsat‐8 spectral imageries, respectively. Four candidate models including stepwise multiple linear regression (SMLR), random forest (RF), cubist (CU) and extreme gradient boosting (XGBoost) Tree were tested in this study. Finally, the digital maps at 30‐m resolution of SOC and SOCD and their uncertainties were prepared using the best‐fit model and the bootstrapping method, respectively. Four soil subgroups (Gypsic Haploxerepts, Typic Calcixerepts, Typic Haploxerepts and Xeric Haplocalcids) were identified across the study area. The covariates had variable contributions on the evaluated models. The XGBoost Tree model generally outperformed other models for prediction of SOC and SOCD (R2 = 0.60, on average). Regardless of soil subgroups, the uncertainty analysis showed that the SOCD map had a low prediction interval range value indicating high accuracy. Additionally, the highest SOCS and SOCD was observed at the top followed by middle and bottom depths in the study area. All subgroups exhibited a decreasing trend of SOCD with increasing depth. A similar trend was also observed for SOCS. The highest SOCD (on average) was observed in Gypsic Haploxerepts (4.71 kg C/m2) followed by Typic Calcixerepts (4.46 kg C/m2), Typic Haploxerepts (4.45 kg C/m2) and Xeric Haplocalcids (4.40 kg C/m2). Overall, the SOCS normalized by area within soil order boundaries was greater in Inceptisols than Aridisols across the study area. The findings of this study provide critical information for sustainable management of soil resources in the area for agricultural production and environmental health in the Miandoab region of Iran.