Volcanic ashfall forecasts are highly dependent on eruption source parameters (ESPs) and synoptic weather conditions at the time and location of the eruption. In New Zealand, MetService and GNS Science have been jointly developing an ashfall forecast system that incorporates four-dimensional high-resolution numerical weather prediction (NWP) and ESPs into the HYSPLIT model, a state-of-the art hybrid Eulerian and Lagrangian dispersion model widely used for volcanic ash. However, these forecasts are based on discrete ESPs combined with a deterministic weather forecast and thus provide no information on output uncertainty. This shortcoming hinders stakeholder decision making, particularly near the geographical margin of forecasted ashfall and in areas with large gradients in forecasted ash deposition. Our study presents a new approach that incorporates uncertainty from both eruptive and meteorological inputs to deliver uncertainty in the model output. To this end, we developed probability density functions (PDFs) for three key ESPs (plume height, mass eruption rate, eruption duration) tailored to New Zealand’s volcanoes and combine them with NWP ensemble datasets to generate probabilistic ashfall forecasts using the HYSPLIT model. We show that the Latin Hypercube Sampling (LHS) technique can be used to representatively span this four-dimensional parameter space and allow us to add uncertainty quantification to rapid response forecast systems. For a case study of a hypothetical eruption at Tongariro, New Zealand we suggest that large parts of New Zealand’s North Island would not receive adequate warning for potential ashfall if uncertainties were not included in the forecasts. We also propose new probabilistic summary products to support public information and emergency responders decision making.