  Numerical solution of the system of partial differential equations arising
from the Reynolds-Averaged Navier-Stokes (RANS) equations with $k-\epsilon$
turbulence model presents several challenges due to the advection dominated
nature of the problem and the presence of highly nonlinear reaction terms.
State-of-the-art software for the numerical solution of the RANS equations
address these challenges by introducing non-differentiable perturbations in the
model to ensure numerical stability. However, this approach leads to
difficulties in the formulation of the higher-order forward/adjoint problems,
which are needed for scalable Hessian-based uncertainty quantification (UQ)
methods. In this note, we present the construction of a UQ-ready flow solver,
i.e., one that is smoothly differentiable and provides not only forward solver
capabilities but also adjoint and higher derivatives capabilities.
