Statistical analysis of the heterogeneity of soil bulk density in regional scale is important to crop growth and tillage management. While most field data focus on bulk density up to a depth of 0.3 m, this study measured the bulk density of 28 in situ soil core locations up to a depth of 1.3 m. Geostatistical method was used to analyze the anisotropic spatial variability and three-dimensional Sequential Gaussian Simulation (SGSIM) was conducted. The results of nugget/sill ratio indicated that the sampling scheme was fine enough to capture the spatial correlation distance of bulk density in the vertical direction. The nugget effect, however, in the horizontal direction suggested that the sampling frequency was not enough for bulk density in that direction. The result of three-dimensional SGSIM is consistent with measured data, indicating that the simulation method is suitable for this research objective. The results may provide basic information for other agricultural management practices.