【Objective】Soil bulk density (ρb) is one of the most important soil physical properties and can be used to characterize soil quality and soil productivity and as a basic parameter to assess soil carbon and nitrogen storage. Furthermore, ρb has a significant effect on transport of soil water and solutes. However, large-scaled soil databases do not encompass much ρb data, owing to the time- and labor-consuming methods to acquire such data, especially the data of deep soil layers in the field. This study was conducted to (1) explore characteristics of regional spatial variation of ρb relative to soil layer across the loess area, (2) determine effects of soil texture, topography, climate and land use on ρb and (3) compare the stepwise regression method with the pedotransfer function method in simulation of spatial variation of ρb.【Method】A total of 243 sampling sites were set based on a grid sampling scheme (40 km×40 km) in the loess area, and ρb of the soils at the depth of 0~10, 10~20 and 20~40 cm of the soil profile and relevant environmental variables, were collected separately, at each sampling site. Spatial variation of ρb was analyzed with the classical statistic method and geostatistical method, respectively. Stepwise regression equation and pedotransfer function equation was used to simulate spatial distribution of soil ρb, separately.【Result】Results show that ρb varied moderately within a soil layer, and generally increased with soil depth in the profile loess area. ρb variability in the 0~40 cm soil layer was moderate according to the coefficient of variation. In general, cropland was the highest in mean ρb, and followed by forestland and grassland. Semivariance of soil ρb of the 0~10, 10~20 and 20~40 cm soil layer can be best fitted by the exponential model, the exponential model and the spherical model, respectively. Soil ρb of the 0~10 cm soil layer exhibited strong spatial dependence and those of the 10~20 and 20~40 cm soil layer did moderate ones. The optimal interval between sampling sites was 5.6~11.2, 70.9~141.7 and 195.4~390.8 km for the 0~10, 10~20 and 20~40 cm soil layer, respectively. Silt content, land use, elevation and slope gradient were the key factors affecting soil ρb in the 0~10 cm layer; silt content, elevation, multi-year mean annual air temperature, aridity and land use were in the 10~20 cm soil layer; and silt content, elevation, land use, multi-year mean annual precipitation, slope gradient and aridity, were in the 20~40 cm soil layer. The pedotransfer function equation explained 38%~52% the variation of ρb, while the stepwise regression equation did only 34%~39%.【Conclusion】Spatial distribution of soil ρb varies significantly with soil depth and vegetation type across the loess area, and is affected jointly by soil texture, topography, climate and land use at the regional scale. The pedostransfer function equation is recommended for modeling and predicting spatial distributions of ρb, particularly for soil layers below 40 cm in the loess area of China.