Synthetic Jet Actuators (SJA) are an Active Flow Control strategy that can decrease the drag caused by boundary layer separation on airplane wings and wind turbine blades. This mechanism consists of a jet that alternately draws in and expels air at a point on the upper surface of the airfoil. The performance of SJAs strongly depends on five control parameters (\textit{i.e.,} jet position and width, amplitude and frequency of oscillations and incidence angle). However, when optimizing such parameters, each of the multiple evaluations of the objective function involves CFD simulations that are computationally very costly. In our approach, we employ a Gaussian Process (GP) surrogate model to substantially reduce the computational burden associated with repeated CFD evaluations. The optimization is conducted using the Efficient Global Optimization (EGO) framework, which iteratively suggests new sampling points by minimizing the Expected Improvement of the objective function(s). In the past, this methodology has been mentioned only briefly for SJA optimization, in the conference presentations [1, 2]. Their work suggested that adding observation noise into the GP can be beneficial to mitigate the uncertainties induced by CFD simulations. In the present study, we provide a more systematic investigation of this approach. We assess the convergence behavior of EGO, and the quality of the resulting optimal designs. The performance of the GP-EGO strategy is compared to our previous work [3], which utilized a Genetic Algorithm and required a substantially large number of CFD evaluations. The comparison is carried out in terms of both the achieved Pareto front and the associated computational cost. Our findings indicate that the surrogate-based method can achieve comparable optimization performance while reducing the number of high-fidelity simulations required.