Currently used assisted history matching algorithms such as differential evolution, particle swarm optimization etc. usually requires a large number of numerical simulation runs in order to converge to acceptable solutions. If each numerical simulation takes a long time to complete, these algorithms become inefficient. In this research, a new assisted history matching tool is presented that can provide multiple solutions of history matching in much less number of numerical simulations. The proposed tool uses Gaussian process based proxy models to provide fast approximate forward solutions which are used in Bayesian optimization to find history match solutions in an iterative manner. The uncertainty in history match solutions is quantified via MCMC sampling on trained forward GP model. In order to check for independence and convergence of the collected MCMC samples, auto-correlation plots and Geweke z-score diagnostic tests are also performed. The converged MCMC samples are then used to quantify uncertainty in EUR of gas reserves via a forecasting GP model. The proposed methodology is successfully applied to a synthetic heterogeneous coalbed methane reservoir. The code and data for this case study is also available online for future studies. The results show that history matching can be performed in approximately four times less number of numerical simulation runs as compared to state of the art differential evolution algorithm. Besides, the P50 estimate of EUR is shown to be in close agreement with truth values for the presented case study.