Abstract: Choosing a suitable prediction model to estimate soil organic matter (SOM) content is not only a prerequisite to improve the accuracy of spatial distribution simulation, but also the basis for monitoring dynamic changes in soil carbon pool and for guiding soil fertility input in farming. In order to achieve this, a research was set up to investigate the advantages of combined traditional Ordinary Kriging (OK) interpolation and Co-Kriging (COK) interpolation in constructing a new model that integrates Cooperative Correlation of auxiliary variables with OK model (CCOK). The following three aspects were thus discussed:1) whether the inclusion of auxiliary variables had an impact on SOM prediction result; 2) what were the differences in SOM prediction results caused by changes in the number of auxiliary information interpolations; and 3) how improved SOM prediction accuracy by cooperative correlation of auxiliary variables. To address these research questions, we collected 329 soil samples from a citrus plantation in Honghuatao Town located in the north Yidu City, Hubei Province. Through physical and chemical analysis, 14 soil properties were extracted. The correlation between SOM and other soil properties were discussed based on Pearson correlation coefficient (r) and available nitrogen was chosen as model auxiliary variable with the most significant correlation with SOM. With reference of OK (the control), we constructed modeling COK (COK₁), global COK (COK₂) and two improved OK models (CCOK₁ and CCOK₂). Among the models, COK₁ was a COK model which used modeling set auxiliary variables to participate in modeling. Based on COK₁, COK₂ changed the modeling set auxiliary variables to global auxiliary variables. CCOK₁ and CCOK₂ represented the OK interpolation models of two forms of functions constructed by the target variables and its auxiliary variables. Some of the results obtained were as follows:1) the range of the nugget/sill proportions of OK, CCOK₁ and CCOK₂ were 25%-75%, which belonged to medium spatial autocorrelation. However, the nugget/sill proportions of COK₁ and COK₂ were less than 25%, belonging to strong spatial autocorrelation. It then showed that the spatial variability of SOM as cross-variance function with auxiliary variables was more easily recognized by semi-variogram models. 2) The predicted SOM in the study area was within 7.38-29.03·kg^-1. Compared with OK interpolation, the strong spatial autocorrelation of COK₁ and COK₂ meant that the spatial distribution of SOM was even more fragmented. Furthermore, plots of CCOK₁ and CCOK₂ predictions were flaky, with digital mapping results of SOM with higher or lower values, which was more consistent with the actual distribution of land use in the study area. 3) The accuracies of COK₁ and OK were similar, but that of COK₂ was higher than the above two. Nevertheless, the correlation coefficients (r) of CCOK₁ and CCOK₂ increased from 0.10 to 0.70 and 0.69, with root mean square errors (RMSE) decreasing by 15.40% and 14.78%, respectively. Finally, the overall accuracy of SOM digital soil mapping was CCOK₁ ≈ CCOK₂ > COK₂ > COK₁ ≈ OK. This indicated that CCOK model minimized error between measured and predicted values in SOM prediction. Thus, the synergy of combined SOM estimation and auxiliary variables was a better correlation than the addition of only auxiliary variables or changing the amount of auxiliary variables. The improved OK model proposed in this study improved the maximum participation of auxiliary information, thereby providing a reliable reference for SOM prediction.