Groundwater numerical models have been widely used as effective tools to analyze and manage water resources. However, the accuracy and reliability of a groundwater numerical model largely depends on model parameters calibration, which is extremely computationally expensive. Therefore, it is highly desirable that efficient optimization algorithms be applied to automatic calibration problems. In this study, we compare the performance of three optimization algorithms and propose a new hybrid method. The algorithms are applied to calibration of a model for part of Beijing water supply. We first outline the three algorithms and briefly describe our hybrid method. The first algorithm referred as PEST in this paper is the Gauss-Marquardt-Levenberg (GML) method including truncated singular value decomposition, which is widely applied in the field of model parameter calibration. As the second one, CMAES_P is a "PEST compatible" implementation of CMA-ES (Covariance Matrix Adaptation Evolution Strategy) global optimization algorithm. PEST derivative-based algorithm and CMAES_P are both encapsulated in the automated parameter optimization software PEST, which has advanced predictive analysis and regularization features to minimize user-specified objective functions. The third one, called Stochastic Radial Basis Function (Stochastic RBF) method, is developed by Regis and Shoemaker (2007), which utilizes radial basis function as the response surface model to approximate the expensive objective function. Our new hybrid method combines Stochastic RBF and PEST derivative-based algorithm, which provides PEST derivative-based algorithm with the starting points found by Stochastic RBF. This paper compares the performances of the aforementioned four algorithms for automatic parameter calibration of a groundwater model on three 28-parameter cases and two synthetic test function calibration problems. We employ the following characteristics as our comparison criteria on all the cases: (1) efficiency in giving good objective function for a given number of function evaluations; (2) performance for different statistical criteria; (3) variability of solutions in multiple trials; (4) improvements if more function evaluations are performed. On the basis of 20 trials, the results indicate that Stochastic RBF is best among the three and CMAES_P is superior to PEST. In addition, our hybrid method still failed to beat Stochastic RBF in highly computationally expensive nonlinear cases. To sum up, our results show that Stochastic RBF method is a more efficient alternative to PEST for automatic parameter calibration of computationally expensive groundwater models. ii