Bulk density (BD) is one of the most important soil physical properties in the hydrological processes of Earth's Critical Zone (CZ). The heterogeneity and related controlling factors of surface BD have received extensive research attention. However, it is unclear whether the effects of controlling factors on BD are homogeneous or vary spatially. This study quantified topsoil BD with a high-density sampling strategy (10 m  10 m; n = 1134) across the ShuangChaGou watershed in the critical zone of the Chinese Loess Plateau (CLP-CZ). Stepwise multiple linear regression (SMLR), regression kriging (RK), random forest (RF) and geographically weighted regression (GWR) were established as modelling approaches to (i) compare the performance of these four models in explaining BD variation, (ii) assess the main factors driving BD, and (iii) evaluate the spatial non-stationary relationship between BD and its explanatory variables. GWR explained the BD variation to the same extent as the RK model (35%), outperforming both the SMLR (15%) and RF models (19%). The variation in topsoil BD depended mainly on: soil organic carbon (SOC) content, normalised difference vegetation index (NDVI), plan curvature, and slope length. The spatial distribution of the local R 2 , regression coefficient magnitude, significance and sign calculated by GWR revealed that the effects of variables accounting for variation in BD varied spatially. The explanatory power of these variables in terms of BD variation was strongest in northern and northeastern areas of the watershed and weakest in southern areas of the watershed. Both significant negative and positive correlations between BD and silt content, NDVI, and plan curvature were found. This study enriched traditional research linked to the complicated relationships between BD and controlling factors by incorporating spatial non-stationarity, which is useful for improving the accuracy of hydrological modelling in the CLP-CZ and in similar CZs worldwide. Highlights• Non-stationary relationships between BD (N = 1134) and its explanatory variables were captured at the watershed scale.