The optimization design of surfactant-enhanced aquifer remediation for dense non-aqueous phase liquids (DNAPLs)-contaminated groundwater is proposed through integrating simulation and optimization models. Many studies have demonstrated that the surrogate model is an effective tool for building a bridge between simulation and optimization models. In this paper, the simulation model was first established to simulate a surfactant-enhanced aquifer remediation process, and on this basis, the ensemble surrogate models were constructed by applying the polynomial regression model, radial basis function neurons network and Kriging surrogate models. Second, the best surrogate model was selected in terms of the three performance indices. Lastly, a non-linear programming optimization model was constructed with the target of minimizing the DNAPLs-contaminated aquifer remediation cost. Meanwhile, the best surrogate model was embedded into the optimization model as a constrained condition, and it was used to reflect the non-linear complex relationship between injection and extraction rates with DNAPLs removal rates instead of simulation models. The results showed that the ensemble surrogate models improved the performance of the single surrogate models. Moreover, ensemble surrogate models improved the computational efficiency, and the optimal strategies have been proved to be an effective guide for contaminant remediation processes.