Wind turbines are critically important in the quest to move towards a non-renewable energy independent world. With the space to add 5M wind turbines, the United States is at the forefront of this transition. Horizontal axis wind turbines (HAWTs) have been studied numerically and experimentally at length. The vast majority of computational fluid dynamics (CFD) studies of HAWTs documented in the open literature have been carried out using two-dimensional simulations. Currently, the available three-dimensional simulations do not provide a comprehensive investigation of the accuracy of different options for modeling of fluid turbulence. In this paper four sets of CFD simulations are carried out using four different turbulence models that are commonly used for engineering level CFD analysis: SST-k-ω, Transition k-kL-ω, Standard k-ε, and Monotonically Integrated Large Eddy Simulation (MILES). These models were compared with experimental performance and coefficient of power results for a small-scale industrial wind turbine with inverse tip speed ratios (λ^(-1) ) in the range 0.072-0.144. They were further investigated to highlight the similarities and differences for the prediction of coefficient of pressure, skin friction coefficient, wall shear, radial velocity and turbulence kinetic energy for λ^(-1) from 0.072 to 0.324. The results showed that no singular model, of the four investigated, was able to consistently predict the power performance with a high degree of accuracy when compared to the experimental results. The models also exhibited both similarities and key differences for the other aspects of flow physics. The results presented in this study highlight the critical role that turbulence modeling plays in the overall accuracy of a CFD simulation, and indicate that end users should be well aware of the uncertainties that arise in CFD results for wind turbine analysis, even when other sources of numerical error have been carefully minimized.