Abstract. Traditional approaches to develop 3D geological models employ a mix of quantitative and qualitative scientific techniques, which do not fully provide quantification of uncertainty in the constructed models and fail to optimally weight geological field observations against constraints from geophysical data. Here, we demonstrate a Bayesian methodology to fuse geological field observations with aeromagnetic and gravity data to build robust 3D models in a 13.5 × 13.5 km region of the Gascoyne Province, Western Australia. Our approach is validated by comparing model results to independently-constrained geological maps and cross-sections produced by the Geological Survey of Western Australia. By fusing geological field data with magnetics and gravity surveys, we show that at 89 % of the modelled region has > 95 % certainty. The boundaries between geological units are characterized by narrow regions with < 95 % certainty, which are typically 400–1000 m wide at the Earth's surface and 500–2000 m wide at depth. Beyond ~ 4 km depth, the model requires drill hole data and geophysical survey data with longer wavelengths (e.g., active seismic) to constrain the deeper subsurface. Our results show that surface geological observations fused with geophysical survey data yield robust 3D geological models with narrow uncertainty regions at the surface and shallow subsurface, which will be especially valuable for mineral exploration and the development of 3D geological models under cover.