A liquid flow microchannel heat sink has been studied and optimized with the help of three-dimensional numerical analysis and multiple surrogate methods. Two objective functions, thermal resistance and pumping power have been selected to assess the performance of the microchannel heat sink. The design variables related to the microchannel top and bottom widths, depth and fin width, which contribute to objective functions, have been identified and design space has been explored through some preliminary calculations. Design of experiments was performed and a three-level full factorial design was selected to exploit the design space. The numerical solutions obtained at these design points were utilized to construct surrogate models namely Response Surface Approximations and Kriging. A hybrid multi-objective evolutionary algorithm coupled with surrogate models and a gradient-based search algorithm is applied to find global Pareto-optimal solutions. Since, the surrogate models are highly problem-dependent, the accuracy of the two surrogate models has been discussed in view of their predictions at on- and off-Pareto-optimal front. The trade-off analysis was performed in view of the two competing objectives. The Pareto-optimal sensitivity (change in value along the Pareto-optimal front) of the design variables has been found out to economically compromise with the design variables contributing relatively less to the objective functions. The application of the multiple surrogate methods not only improves quality of multi-objective optimization but also gives the feedback of the fidelity of the model near the optimum region.