A first study on the continuous adjoint formulation for aerodynamic optimization design of high pressure turbines based on S2 surface governed by the Euler equations with source terms is presented. The objective function is defined as an integral function along the boundaries, and the adjoint equations and the boundary conditions are derived by introducing the adjoint variable vectors. The gradient expression of the objective function then includes only the terms related to physical shape variations. The numerical solution of the adjoint equation is conducted by a finite-difference method with the Jameson spatial scheme employing the first and the third order dissipative fluxes. A gradient-based aerodynamic optimization system is established by integrating the blade stagger angles, the stacking lines and the passage perturbation parameterization with the quasi-Newton method of Broyden–Fletcher–Goldfarb–Shanno (BFGS). The application of the continuous adjoint method is validated through a single stage high pressure turbine optimization case. The adiabatic efficiency increases from 0.8875 to 0.8931, whilst the mass flow rate and the pressure ratio remain almost unchanged. The optimization design is shown to reduce the passage vortex loss as well as the mixing loss due to the cooling air injection.