A design method was developed for automated, systematic design of hydrokinetic turbine rotor blades. The method coupled a Computational Fluid Dynamics (CFD) solver to estimate the power output of a given turbine with a surrogate-based constrained optimization method. This allowed the characterization of the design space while minimizing the number of analyzed blade geometries and the associated computational effort. An initial blade geometry developed using a lifting line optimization method was selected as the base geometry to generate a turbine blade family by multiplying a series of geometric parameters with corresponding linear functions. A performance database was constructed for the turbine blade family with the CFD solver and used to build the surrogate function. The linear functions were then incorporated into a constrained nonlinear optimization algorithm to solve for the blade geometry with the highest efficiency. A constraint on the minimum pressure on the blade could be set to prevent cavitation inception.