Groundwater contamination source identification (GCSI) is critical for taking effective measures to protect groundwater resources, assess risks, mitigate disasters, and design remediation strategies. Simulation-optimization techniques have been effective tools for GCSI. However, previous studies have applied individual surrogate models when replacing simulation models, rather than making efforts to combine various methods to improve the approximation accuracy of the surrogate model over the simulation model. In this study, the kernel extreme learning machine (KELM) model was proposed to enhance the surrogate model, and to approach GCSI problems, especially those of dense nonaqueous phase liquid-contaminated aquifers, more effectively. In addition, a kriging model and a support vector regression (SVR) model were built and compared with the KELM model, and various ensemble surrogate (ES) modeling techniques were applied to establish four ES models. Results showed that the KELM model was more accurate than the kriging and SVR models; however, the ES models performed much better than the three individual surrogate models. The most precise ES model increased the certainty coefficient (R 2 ) to 0.9837, whereas limiting the maximum relative error to 13.14%. Finally, a mixedinteger nonlinear programming optimization model was established to identify the groundwater contamination source in terms of location and release history, and simultaneously assess aquifer parameters. The ES model developed in this article could reasonably predict the system response under given operation conditions. Replacement of the simulation model by the ES model considerably reduced the computation burden of the simulation-optimization process and simultaneously achieved high computation accuracy.