In this paper, we present an approach to inversely quantify the uncertainty of MCFD simulations through a data-driven modular Bayesian inference. Both the model parameter uncertainty and the model form uncertainty are evaluated in the proposed approach. Considering the high-dimensionality of parameter space related to the solver, we performed a sensitivity analysis to reduce the input parameter dimension for the Bayesian inference. Based on the reduced parameter dimension, surrogate models based on Gaussian Process (GP) and Principal Component Analysis (PCA) are constructed to reduce the computational cost in the Bayesian inference. Two case studies based on the proposed approach are performed, focusing on two-phase flow dynamics and on wall boiling heat transfer, respectively. In case study I, we are able to construct a GP-based surrogate model based on 8 principal components to represent the total 208 MCFD solver outputs. Moreover, both cases show that the proposed approach is able to quantify and reduce the parameter uncertainties with the support of experimental measurements. The posterior uncertainties of investigated parameters have 50%-90% narrowed uncertainty ranges compared to their prior uncertainties. Furthermore, a forward uncertainty propagation of the MCFD solver with the obtained uncertainties shows that the agreement between the solver predictions and experimental measurements are significantly improved. (C) 2019 Elsevier Ltd. All rights reserved.