In this thesis, detailed uncertainty quantification studies focusing on the closure coefficients of eddy-viscosity turbulence models for several flows using two CFD solvers have been performed. Three eddy viscosity turbulence models considered are: the one-equation Spalart-Allmaras (SA) model, the two-equation Shear Stress Transport (SST) k-[omega] model, and the one-equation Wray-Agarwal (WA) model. OpenFOAM and ANSYS Fluent are used as flow solvers. Uncertainty quantification analyses are performed for subsonic flow over a flat plate, subsonic flow over a backward-facing step, and transonic flow past an axisymmetric bump. In the case of flat plate, coefficients of pressure, lift, drag, and skin friction are considered to be the output quantities of interest. In case of the backward-facing step, these quantities are considered along with the separation bubble size. In case of an axisymmetric transonic bump, the drag coefficient, lift coefficient, separation point and reattachment point are considered. In addition to these four quantities, global uncertainty is employed on every node in the flow for Reynolds shear stress to determine which areas of the flow the closure coefficients contribute most to the uncertainty. Uncertainty quantification is conducted using DAKOTA developed by Sandia National Laboratories using stochastic expansions based on non-intrusive polynomial chaos. All closure xii coefficients are treated as epistemic uncertain variables, each defined by a specified range. The influence of the closure coefficients on output quantities is assessed using the global sensitivity analysis based on variance decomposition. This yields Sobol indices which are used to rank the contributions of each constant. A comparison of the Sobol indices between the turbulence models, flow cases, and flow solvers is conducted. This research identifies closure coefficients for each turbulence model that contribute significantly to uncertainty in the model predictions; this information can then be used to improve the prediction capability of the models in separated flow region by a more judicious choice of the closure coefficients.