Given spatial variability of soil organic carbon (SOC) as a key indicator of soil quality, its precise regional mapping is essential. This study compared random forest (RF), regression kriging (RK), and extreme gradient boosting (XGBoost) algorithms to identify the most influential factors controlling SOC distribution in the Karun-3 watershed. 158 composite surface soil samples (0–5 cm) were collected using a conditioned Latin Hypercube Sampling (cLHS) design. Standard laboratory methods measured SOC, sand, silt, clay, electrical conductivity (EC), pH, calcium equivalent carbonate (CCE), and bulk density (BD). Additionally, 31 auxiliary variables—topographic, spectral, climatic indices, and soil characteristics—were used to model SOC spatial variability. These variables were pre-processed by removing those with near-zero variance and those highly intercorrelated (r ≥ 0.8) with a variance inflated factor (VIF) > 10. The XGBoost algorithm outperformed the others, demonstrating the highest coefficient of determination (R2 = 0.628), the lowest error metrics (RMSE = 0.473 and MAE = 0.316), and the highest ratio of performance to interquartile range (RPIQ = 1.98). BD, CCE, Carbonate Index (Carb.Ind), Normalized Difference Moisture Index (NDMI), and Length-Slope Factor (LS-Factor) contributed most significantly to SOC distribution. SOC distribution maps indicated that the northern and northwestern parts of the area contain the highest levels, while the central, southern, and southeastern parts have the lowest. Uncertainty maps revealed that predictive map values are less reliable in parts of the northwestern, eastern, and southern study area. Overall, simultaneous application of predictive and uncertainty maps is necessary to understand both predicted values and their reliability.