With recent advances in high-performance computing, computational fluid dynamics (CFD) modeling has become an integral part in the engineering analysis and even in the design process of marine vessels and propulsors. In aircraft wing design, CFD has been integrated with numerical optimization and adjoint methods to enable high-fidelity aerodynamic shape optimization with respect to large numbers of design variables. There is a potential to use some of these techniques for maritime applications, but there are new challenges that need to be addressed to realize that potential. This work presents a solution to some of those challenges by developing a CFD-based hydrodynamic shape optimization tool that considers cavitation and a wide range of operating conditions. A previously developed three-dimensional compressible Reynolds-averaged Navier‐Stokes (RANS) solver is extended to solve for nearly incompressible flows, using a low-speed preconditioner. An efficient gradient-based optimizer and the adjoint method are used to carry out the optimization. The modified CFD solver is validated and verified for a tapered NACA 0009 hydrofoil. The need for a large number of design variables is demonstrated by comparing the optimized solution obtained using different number of shape design variables. The results showed that at least 200 design variables are needed to get a converged optimal solution for the hydrofoil considered. The need for a high-fidelity hydrodynamic optimization tool is also demonstrated by comparing RANS-based optimization with Euler-based optimization. The results show that at high lift coefficient (C L) values, the Euler-based optimization leads to a geometry that cannot meet the required lift at the same angle of attack as the original foil due to inability of the Euler solver to predict viscous effects. Single-point optimization studies are conducted for various target C L values and compared with the geometry and performance of the original NACA 0009 hydrofoil, as well as with the results from a multipoint optimization study. A total of 210 design variables are used in the optimization studies. The optimized foil is found to have a much lower negative suction peak, and hence delayed cavitation inception, in addition to higher efficiency, compared to the original foil at the design C L value. The results show significantly different optimal geometry for each C L, which means an active morphing capability was needed to achieve the best possible performance for all conditions. For the single-point optimization, using the highest C L as the design point, the optimized foil yielded the best performance at the design point, but the performance degraded at the off-design C L points compared to the multipoint design. In particular, the foil optimized for the highest C L showed inferior performance even compared to the original foil at the lowest C L condition. On the other hand, the multipoint optimized hydrofoil was found to perform better than the original NACA 0009 hydrofoil over the entire operation profile, where the overall efficiency weighted by the probability of operation at each C L, is improved by 14.4%. For the multipoint optimized foil, the geometry remains fixed throughout the operation profile and the overall efficiency was only 1.5% lower than the hypothetical actively morphed foil with the optimal geometry at each C L. The new methodology presented herein has the potential to improve the design of hydrodynamic lifting surfaces such as propulsors, hydrofoils, and hulls.