Using turbulence models with finite element methods (FEM) can be challenging as the turbulence variables can assume negative non-physical values and hinder solution convergence. A modified k–ω model was recently proposed by Stefanski et al. (2018) to be used with finite element solvers of compressible flows. The model overcomes this issue by replacing k and ω with working variables that ensure positivity and smoothness of k and ω. In this work the applicability of this model for high-order FEM simulations of incompressible flows was examined. The model was implemented for incompressible flow in an hp-FEM solver using streamline Petrov-Galerkin discretization and was validated and verified using a fully-developed channel flow and a boundary layer flow over a flat plate. Several aspects of the turbulence model behavior were studied. These included the possibilitty of getting orders of accuracy higher than 2, and the model’s sensitivity to freestream values of k and ω. The results suggested that higher orders of accuracy are possible when quadratic and quartic basis functions are used. However, this depended on the way the boundary condition for ω was defined. The commonly used boundary condition for ω, which depends on the wall-distance of the first grid point resulted in poor orders of accuracy compared to the so-called slightly-rough-surface boundary condition which is independent of the wall distance of the first grid point. Additionally, results indicated that increasing the nondimensional wall distance of the first gridpoint makes it more sensitive to the value of ω on the wall. Adding a cross-diffusion term to the transport equation for ω is known to significantly improve the accuracy of turbulence model prediction for certain flows and reduce the sensitivity of the original k–ω model to freestream values of turbulence variables. Following a more recent version of k–ω model, this term was added to the turbulence model and some other modifications including a different production term with a stress-limiter were applied. The drag coefficient of the flat plate from the new turbulence model showed similar sensitivity to the freestream values of turbulence variables as the model of Stefanski et al. (2018).