The creation of a contour map of the water table in an unconfined aquifer based on head measurements is often the first step in any hydrogeological study. Geostatistical interpolation methods (e.g., kriging) may provide exact interpolated groundwater levels at the measurement locations but often fail to represent the hydrogeological flow system. A physically based, numerical groundwater model with spatially variable parameters and inputs is more adequate in representing a flow system. Because of the difficulty in parameterization and solving the inverse problem, however, a considerable difference between calculated and observed heads will often remain. In this study the water‐table interpolation methodology presented by Fasbender et al. (2008), in which the results of a kriging interpolation are combined with information from a drainage network and a digital elevation model (DEM), using the Bayesian data fusion framework, is extended to incorporate information from a tuned analytic element groundwater model. The resulting interpolation is exact at the measurement locations whereas the shape of the head contours is in accordance with the conceptual information incorporated in the groundwater‐flow model. The Bayesian data fusion methodology is applied to a regional, unconfined aquifer in central Belgium. A cross‐validation procedure shows that the predictive capability of the interpolation at unmeasured locations benefits from the Bayesian data fusion of the three data sources (kriging, DEM, and groundwater model), compared to the individual data sources or any combination of two data sources.