Probabilistic volcanic hazard assessments require (1) an identification of the hazardous volcanic source; (2) estimation of the magnitude-frequency relationship for the volcanic process; (3) quantification of the dependence of hazard intensity on magnitude and external conditions; and (4) estimation of hazard exceedance from the magnitude-frequency and hazard intensity relationship. For volcanic mass flows, quantification of the hazard intensity is typically undertaken through the use of computationally expensive mass flow simulators. However, this computational expense restricts the number of samples that can be used to produce a probabilistic assessment and limits the ability to rapidly update hazard assessments in response to (e.g.) changing source probabilities. We develop an alternate approach to defining hazard intensity through a surrogate model that provides a continuous estimate of simulation outputs at negligible computational expense, demonstrated through a probabilistic hazard assessment of dome collapse (block-and-ash) flows at Taranaki volcano, New Zealand. A Gaussian Process emulator trained on a database of simulations is used as the surrogate model of hazard intensity across the input space of possible dome collapse volumes and configurations, which is then sampled using a volume-frequency relationship of dome collapse flows. The demonstrated technique is a tractable solution to the problem of probabilistic volcanic hazard assessment, with the surrogates providing a good approximation of the simulator at very limited computational expense, and is generally applicable to volcanic hazard and geo-hazard assessments that are limited by the demands of numerical simulations.