Summary: An adaptive finite element method for solving incompressible turbulent flows using the \(k\)-\(\varepsilon\) model of turbulence is presented. Solutions are obtained in primitive variables using a highly accurate quadratic finite element on unstructured grids. Two error estimators are presented that take into account in a rigorous way the relative importance of the errors in velocity, pressure, turbulence variables, and eddy viscosity. The efficiency and convergence rate of the methodology are evaluated by solving problems with known analytical solutions. The method is then applied to turbulent free shear flows, and predictions are compared to measurements.