Objective The aims were to precisely map the spatial distribution of soil organic carbon (SOC), in order to widely carry out the precision fertilization for sustainable agricultural development.

Method A zoning regression kriging (ZRK) method considered spatial heterogeneity and interactive environmental variables was proposed to draw a high-precision spatial distribution map of SOC content in northern France. The global ordinary kriging (OK), co-kriging (COK), regression kriging (RK), random forest (RF), and OK, COK, and RK based on one single environmental variable and geographical detector were used for comparison.

Result ① The SOC content in the study area was driven by the interaction of various environmental variables, and the driving effects of different spatial homogeneous zones were different. ② The spatial hierarchical clustering method could effectively use the spatial homogeneous zones of SOC content. Compared with the global OK, COK and RK models, the R² values of the ZOK, ZCOK and ZRK methods based on spatial hierarchical clustering partitions were improved by 32% ~ 36 %. ③ The kriging methods based on single environmental variable and geographical detector partition performed generally in this study, and did not effectively improve the SOC prediction accuracy. ④ RF performed poorly in the study, with lower R² and higher RMSE and MAE, while ZRK performed best in all models, with the highest R² and the lowest RMSE and MAE.

Conclusion The ZRK method, which takes into account spatial heterogeneity and interactive environmental variables, can effectively use the spatial homogeneous zones of SOC content and draw the spatial distribution map of SOC content considered the interaction of environmental variables, which provides a new model reference for digital soil mapping in spatial heterogeneous regions.