Abstract: Abstract: The spatial variability of soil organic carbon (SOC) is one of the reasons leading to uncertainty in the estimation of carbon stocks. Simulation and analysis research of SOC spatial distribution, especially the three-dimensional (3D) spatial distribution characteristics, is of great significance for revealing the soil nutrient and pollutant migration, the precise management of farmland and agricultural sustainable development. The 3D spatial distribution pattern of SOC was estimated using 3D ordinary Kriging with regard to the isotropy (3DOKI), 3D ordinary Kriging and 3D CoKriging coupling Markov with regard to anisotropy (3DOKA and 3DCKMA) in the region. Based on estimation results, the paper revealed 3D distribution pattern of regional SOC content. The mean root mean square error (RMSE), relative improvement (RI), the scatter diagram of SOC simulated and measured value and specific value coverage ratio were used to estimate the accuracy of different approaches, and the mean squared deviation ratio (MSDR) and the 3D variogram were used to evaluate model fitting effect and spatial local uncertainty. The results showed that 3D spatial distribution patterns of SOC content for the 3 kinds of estimation methods were basically the same, and the overall spatial distribution pattern was SOC content of the west and the north was higher, and that of the east and the south was relatively low. SOC contents of 5 layers (0-20, >20-40, >40-60, >60-80 and >80-100 cm) were 11.88±5.76, 10.08±4.89, 8.40±5.49, 7.83±5.89 and 7.17±5.22 g/kg, respectively. With the increase of soil depth, the SOC content of the patch gradually reduced. No matter which kind of method, the spatial distribution of SOC in different soil depth was similar, that was, the distribution characteristics of surface layer of the patch were also the embodiment of the deep soil layer. Compared to 3DOKI, the spatial search strategy and its parameters with regard to anisotropy were able to reduce defect (Bovine and central tendency effect) of Kriging method in a certain extent. The RMSE values for 3DOKI, 3DOKA and 3DCKMA were 3.1645, 2.0523 and 1.6215, respectively, so the RMSE value was the minimum for 3DCKAM, whose RI reached nearly 50%. Specific value coverage ratio for 3DOKI, 3DOKA and 3DCKMA were 33.12%, 57.83% and 76.15%, respectively, which was the largest for 3DCKMA method, whose MSDR value (1.4409) was the most close to 1, and variance was the least. The same position CoKriging method coupling Markov was more accurate, and its model fitting effect was the best, and uncertainty was the smallest, which could better highlight the volatility and reflect the specific value. The related research results of this paper will provide method reference for the regional soil properties, and provide technical guidance for the reasonable and scientific research of the spatial distribution pattern of regional SOC.