Numerical shape optimization is a key approach for improving the aerodynamic performance of turbomachinery components, enabling gains in efficiency, stability margin, and loss reduction. These optimization problems rely on the steady RANS equations closed by a turbulence model and involve many geometric parameters, making the computation of shape gradients a central issue for gradient-based optimization. In industrial aerodynamic applications, the number of shape parameters largely exceeds the number of objective functionals, resulting in a computationally prohibitive cost for direct differentiation. The adjoint method overcomes this limitation by enabling the computation of all shape gradients at a cost independent of the number of parameters. This advantage comes at the expense of solving a large, stiff and ill-conditioned adjoint linear system resulting from the full differentiation of the discrete RANS solver, including the turbulence model for gradient accuracy improvement [1][2]. Within the elsA solver [3], adjoint-based shape optimization has so far been limited to the Spalart–Allmaras turbulence model, a popular choice for turbomachinery shape [4][5]. While its single-equation formulation provides simplicity and robustness, it shows limited accuracy for complex internal turbomachinery flows [6]. To address this limitation, the present work focuses on the two-equation Smith k-l turbulence model [7], which offers a more detailed representation of turbulent mechanisms while remaining suitable for internal flows. This paper presents, for the first time, the complete differentiation of the Smith k-l turbulence model within the steady RANS adjoint solver of elsA. The resulting shape gradients are validated against finite-difference approximations on a two-dimensional backward-facing step. A shape optimization is then performed on this configuration using the k-l model, and the results are compared with those obtained using the Spalart–Allmaras model. This comparison highlights the benefits and limitations of the k-l model for adjoint-based shape optimization of internal flows.