Interest in simulation-based optimal design methods is growing rapidly in industry, as companies need to maximize performance and efficiency while minimizing development costs and development times. This is particularly true for transient applications characterized by complex, violent free surfaces, such as Pelton turbines. Conventional mesh-based methods struggle with these applications due to the extensive mesh adaptation requirements and the challenges of adequately tracking the evolving free surfaces. In this context, particle-based, mesh-free Lagrangian numerical methods, such as SPH, offer significant advantages. Although simulations using SPH methods are already fully integrated into industrial workflows, the use of an adjoint SPH method for shape optimization is still unknown. This work presents a novel numerical method for fluid-dynamic applications of unsteady adjoint shape optimization methods computed with SPH. To close this gap, we describe the derivation and implementation of a weakly compressible, unsteady, primal[1]/continuous-adjoint Riemann SPH solver that can be combined with an Arbitrary Lagrangian-Eulerian formulation of the governing equations (SPH-ALE) to determine the design sensitivities of the cost functional to be minimized. The corresponding adjoint equations are derived from the weakly compressible isentropic Euler equations, discretized using the Riemann SPH-ALE method. The latter employs an efficient reduced-order model to provide the primal flow during the time-reversed adjoint computation [2]. Examples included refer to generic and practical applications using a force-based cost functional.