Over the last two decades, the focus of research in multiphase CFD has been on the development of advanced numerical methods, which leveraged the rapid increase in computing power. Even with the advent of exascale computing in the foreseeable future, detailed simulations of industrial-scale multiphase flows will remain out of reach for decades to come. Within many engineering processes, such as fluidized bed reactors, two-phase flow instabilities often lead to ‘demixing’ resulting in spatially non-uniform suspensions that obstruct chemical/thermal efficiencies, which current turbulence models fail to capture at large scales. As highly-resolved multiphase flow data continues to come online, new techniques are needed to integrate this information across scales. In addition, to aid in decision making, multiphase flow simulations have to be augmented to estimate underlying uncertainties in simulation components. In this talk, we will present a data-driven framework for model closure of the multiphase Reynolds Average Navier—Stokes (RANS) equations. Data generated from high-fidelity simulations are used in combination with state-of-the-art inverse modeling and machine learning techniques to (i) quantify model form uncertainty in existing models and (ii) infer the functional form of new turbulence models across a broad range of two-phase flow regimes.