Abstract: To address the high computational cost of high-fidelity numerical simulations used in the aerodynamic optimization of vertical-axis wind turbine (VAWT) airfoils, an optimization framework integrating XFOIL, Gaussian process regression (GPR), and a genetic algorithm (GA) was developed using the NACA 0021 airfoil as the baseline. The leading-edge radius ratio, maximum thickness ratio, maximum camber, and trailing-edge flap deflection angle were selected as design variables. A total of 480 parametric airfoils were generated, and their aerodynamic performance was evaluated using XFOIL at Ma = 0.03, Re = 1.2 × 10⁵, and α = 8° to establish a lift-to-drag ratio dataset. Of these samples, 380 were used for training and 100 for testing. Four GPR models employing squared exponential (SE), rational quadratic (RQ), Matérn 3/2, and Matérn 5/2 kernels were constructed, and their kernel hyperparameters were optimized using GA. The results showed that the RQ kernel achieved the best overall predictive performance. Following optimization, the root mean square error (RMSE) decreased from 3.046 to 2.571, while the coefficient of determination (R²) increased from 0.944 to 0.956. The Optimized 0021 airfoil was subsequently obtained by applying GA to identify the geometric parameter combination that maximized the lift-to-drag ratio. CFD validation showed that, at α = 8°, the lift-to-drag ratio of the optimized airfoil was 6.4% higher than that of the baseline airfoil, and the relative error between the GA-GPR prediction and the CFD result was 4.4%. Over the angle-of-attack range of 0°–20°, the maximum lift coefficient and maximum lift-to-drag ratio increased by 14.6% and 12.3%, respectively. The proposed method achieves a favorable balance between predictive accuracy and computational efficiency and provides a reference for the aerodynamic optimization of VAWT airfoils.