Calibration of a multi-phase, geothermal flow model is a crucial step for matching model outputs to available data and hence selecting a good set of parameters for a certain model. In general, it requires a large number of evaluations of an expensive computer simulator. Due to the high cost of each flow model evaluation it is essential that the optimization method selected is able to converge quickly to a good solution. It is well known that nonlinear models may present several local minima, and hence it is possible to get stuck in one of these low-quality solutions using derivative-based optimization techniques only. Global optimization methods are not very efficient, mainly due to the extremely high number of evaluations required for convergence. In this paper we explore the use of a global, derivative-free, optimization method linked with a surrogate model for the calibration of a geothermal reservoir model. The use of a surrogate model allows us to explore the input space in a faster way than using the high fidelity flow simulator. The surrogate model is subsequently trained with new samples suggested from the global method. This workflow is tested in an actual geothermal reservoir model which is complex enough to simulate conditions observed in real life problems. The use of this approach allowed us to calibrate a highly parameterized numerical model using a robust methodology which could be applied to other similar subsurface flow problems.