History matching is commonly performed in reservoir simulations to calibrate model parameters and to improve prediction accuracy. History matching problems often have non-unique solutions, i.e., there exist different combinations of parameter values that all yield the simulation results matching the measurements. In such a situation, finding a single solution matching the observations does not guarantee a correct prediction for future production. Alternatively, a more reliable prediction should be made with an uncertainty quantification based on all different possible scenarios of the model parameters. Bayesian theorem provides a theoretical foundation to represent different solutions and to quantify the uncertainty with the posterior probability density function (PDF). Lacking an analytical expression, the posterior PDF are often shown with a sample of realizations, each representing a possible scenario. This paper presents a novel sampling algorithm aiming to deal with two commonly encountered difficulties in the sampling process. First, a typical sampling method requires intensive model evaluations and hence may cause unaffordable computational burden. To alleviate this burden, our algorithm uses a Gaussian process (GP)-based surrogate as an approximation of the computationally expensive reservoir model to speed up the sampling process. The GP surrogate is adaptively refined locally such that the necessary approximation accuracy is achieved with a minimum level of computational cost. Secondly, when the dependent relationship between observation variables and input parameters is nonlinear, the posterior PDF could be in a complex form, such as multimodal, which is difficult to sample from. To tackle this difficulty, a Gaussian mixture model (GMM) is used as the proposal PDF to explore the parameter space. The GMM is flexible to approximate different distributions and is particularly efficient when the posterior is multimodal. The developed approach is tested with an illustrative history matching problem and shows its capability in handling the above-mentioned issues. Multimodal posterior of the testing problem is captured and are used to give a reliable production prediction with uncertainty quantification. The new algorithm reveals a great improvement in terms of computational efficiency comparing previously studied approaches for the sample problem.