As all kinds of physics-based and data-driven models are emerging in the fields of hydrologic and hydraulic engineering, Bayesian model averaging (BMA) is one of the popular multi-model methods used to account for the various uncertainty sources in the flood modeling process and generate robust ensemble predictions based on multiple competitive candidate models. The reliability of BMA parameters (weights and variances) determines the accuracy of BMA predictions. However, the uncertainty in the BMA parameters with fixed values, which are usually obtained from the Expectation-Maximization (EM) algorithm, has not been adequately investigated in BMA-related applications over the past few decades. Given the limitations of the commonly used EM algorithm, the Metropolis-Hastings (M-H) algorithm, which is one of the most widely used algorithms in the Markov Chain Monte Carlo (MCMC) method, is proposed to estimate the BMA parameters and quantify their associated uncertainty. Both numerical experiments and the one-dimensional HEC-RAS models are employed to examine the applicability of the M-H algorithm with multiple independent Markov chains. The performances of the EM and M-H algorithms in the BMA analysis are compared based on the daily water stage predictions from 10 model configurations. The results show that the BMA weights estimated from both algorithms are comparable, while the BMA variances obtained from the M-H MCMC algorithm are closer to the given variances in the numerical experiment. Moreover, the normal proposal distribution used in the M-H algorithm can yield narrower distributions for the BMA weights than those from the uniform prior. Overall, the MCMC approach with multiple chains can provide more information associated with the uncertainty of BMA parameters and its prediction performance is better than the default EM algorithm in terms of multiple evaluation metrics as well as algorithm flexibility.