Lava flow models have been created to forecast the behavior of lava in order to assess their hazard to communities near volcanoes. Many flow model frameworks are deterministic, where flow parameters are assumed to be well known and model output produces a single estimated solution of the flow characteristic in question (e.g. runout length, areal extent of inundation). Probabilistic approaches to modeling lava flows may have inherent advantages in flow forecasting as these approaches consider parameter uncertainty and provide model results that are not binary. For instance, a deterministic lava flow model might provide a "safe" or "not safe" result for any location on a map. A probabilistic flow model might instead provide a range of likelihoods of inundation given a lava flow; locations will have a 0- 100% chance of being engulfed in lava. Here, a new Cellular Automata lava flow program named MOLASSES is used to generate probabilistic hazard maps. MOLASSES simulates the ultimate areal extent of a lava flow deterministically using source location(s), a lava flow thickness estimate, a total flow volume, and a digital elevation model. To create a probabilistic model of flow inundation, each of these parameters are varied based on their uncertainties or on statistical models from previous flows. The simulator is then run as a Monte Carlo, which is possible due to the speed of the code. Two types of probabilistic hazard maps are created. First, inundation probability for the 2012-3 lava flows at Tolbachik, Kamchatka, Russia, is assessed and compared to the real lava flow. Second, a long-term flow inundation forecast is created for the East Snake River Plain Volcanic Field in Idaho, USA, using a database of known locations and morphologies of past flows. The resulting probabilistic maps can potentially be used in hazards assessment both long-term and during a crisis.