Axial waterjets are widely used for marine propulsion due to their efficiency and maneuverability. However, conventional design procedures heavily rely on empirical correlations and simplified models, limiting their ability to fully exploit the hydrodynamic performance potential of these devices. The study highlights how Simulation-Based Design Optimization (SBDO) approaches, coupled with the high-fidelity simulations required to hydrodynamically characterize the complex phenomena that occur in the case of waterjets, can enable the identification of non-intuitive design improvements over a wider design space that may be missed by traditional methods. In particular, the Reynolds-Averaged Navier-Stokes (RANS) equations are used to provide accurate performance predictions, capturing complex flow phenomena such as secondary flows (i.e., leakage vortices) and pressure distributions critical to waterjet design, of systematically varied configurations using a 42-dimensional parametric model. Simplified key performance indicators, in the specific cavitation inception obtained from the non-cavitating analysis, work in conjunction with the calculated hydraulic efficiency to identify geometries capable of improving (or not worsening) efficiency while postponing cavitation. The systematic and automated analysis of thousands of different configurations, iteratively modified by a genetic algorithm, is finally able to identify better waterjets, whose performances are confirmed by dedicated cavitating RANSE analyses. This demonstrates how RANS-based simulations, integrated with optimization algorithms, can lead to superior axial waterjet designs, providing a flexible, more robust, and effective methodology compared to conventional approaches.