To provide soil hydraulic parameters for simulating soil hydrological processes, classical statistics were used to characterize the spatial distribution of soil bulk density (BD) within different soil depths (0~10cm, 10~20cm and 20~40cm) along south to north transect of the Loess Plateau. Furthermore, the multiple stepwise regressions, pedotransfer functions and the first order autoregressive statespace models were applied to simulate the distribution of BD. The magnitude of BD variability in the 0~20cm soil layer was moderate according to the coefficient of variations, while the BD variability in the 20~40cm soil layer was weak. The key factors affecting the spatial distribution of BD differed in different soil depths. Soil organic carbon, clay and sand contents were the key factors in the 0~10cm soil layer;soil organic carbon content, clay and sand contents and precipitation were the key factors in the 10~20cm soil layer;while in the 20~40cm soil layer, clay and sand contents, precipitation and land use were the key factors to affect the spatial distribution of BD along south to north transect. State-space models were consistently more effective than multiple stepwise regression functions and pedotransfer functions for estimating spatial distribution of BD. State-space model that included clay and sand contents, precipitation and land use factors showed the best simulation result, and the combination of such variables explained 92.3% of the total variation of BD. State-space models were recommended for studying spatial relations between soil bulk density and other variables on the Loess Plateau.