Parameterised models that predict the gravitational-wave (GW) signal from merging black holes are used to extract source properties from GW observations. The majority of research in this area has focused on developing methods capable of producing highly accurate, point-estimate, predictions for the GW signal. A key element missing from every model used in the analysis of GW data is an estimate for how confident the model is in its prediction. This omission increases the risk of biased parameter estimation of source properties. Current strategies include running analyses with multiple models to measure systematic bias however, this fails to accurately reflect the true uncertainty in the models. In this work we develop a probabilistic extension to the phenomenological modelling workflow for non-spinning black holes and demonstrate that the model not only produces accurate point-estimates for the GW signal but can be used to provide well-calibrated local estimates for its uncertainty. Our analysis highlights that there is a lack of Numerical Relativity (NR) simulations available at multiple resolutions which can be used to estimate their numerical error and implore the NR community to continue to improve their estimates for the error in NR solutions published. Waveform models that are not only accurate in their point-estimate predictions but also in their error estimates are a potential way to mitigate bias in GW parameter estimation of compact binaries due to unconfident waveform model extrapolations.