  The Reynolds-averaged Navier-Stokes (RANS) equations for steady-state assessment of incompressible turbulent flows remain the workhorse for practical computational fluid dynamics (CFD) applications. Consequently, improvements in speed or accuracy have the potential to affect a diverse range of applications. We introduce a machine learning framework for the {surrogate modeling of steady-state turbulent eddy viscosities for RANS simulations, given the initial conditions. This modeling strategy} is assessed for parametric interpolation, while numerically solving for the pressure and velocity equations to steady state, thus representing a framework that is hybridized with machine learning. We achieve {competitive} steady-state results with a significant reduction in solution time when compared to those obtained by the Spalart-Allmaras one-equation model. This is because the proposed methodology allows for considerably larger relaxation factors for the steady-state velocity and pressure solvers. Our assessments are made for a backward-facing step with considerable mesh anisotropy and separation to represent a practical CFD application. For test experiments with \textcolor{black}{either} varying inlet velocity conditions or step heights we see time-to-solution reductions around a factor of 5. The results represent an opportunity for the rapid exploration of parameter spaces that prove prohibitive when utilizing turbulence closure models with multiple coupled partial differential equations. \blfootnote{Code available publicly at \texttt{https://github.com/argonne-lcf/TensorFlowFoam}}. 
A turbulent eddy-viscosity surrogate modeling framework for   Reynolds-Averaged Navier-Stokes simulations