Numerical modeling is a significant tool to understand the dynamic characteristics of contaminants transport in groundwater. The automatic calibration of highly parametrized and computationally intensive numerical models for the simulation of contaminant transport in the groundwater flow system is a challenging task. While existing methods use general optimization techniques to achieve automatic calibration, the large numbers of numerical model evaluations required in the calibration process lead to high computing overhead and limit the efficiency of model calibration. This paper presents a Bayesian optimization (BO) method for efficient calibration of numerical models of groundwater contaminant transport. A Bayes model is built to fully represent calibration criteria and derive the objective function for model calibration. The efficiency of model calibration is made possible by the probabilistic surrogate model and the expected improvement acquisition function in BO. The probabilistic surrogate model approximates the computationally expensive objective function with a closed-form expression that can be computed efficiently, while the expected improvement acquisition function proposes the most promising model parameters to improve the fitness to the calibration criteria and reduce the uncertainty of the surrogate model. These schemes allow us to find the optimized model parameters effectively by using a small number of numerical model evaluations. Two case studies for the calibration of the Cr(VI) transport model demonstrate that the BO method is effective and efficient in the inversion of hypothetical model parameters, the minimization of the objective function, and the adaptation of different model calibration criteria. Specifically, this promising performance is achieved within 200 numerical model evaluations, which substantially reduces the computing budget for model calibration.