This paper describes the development of a numerical model for studying the evolution of a wave train, shoaling and breaking in the surf zone. The model solves the Reynolds equations for the mean (ensemble average) flow field and the k–ε equations for the turbulent kinetic energy, k, and the turbulence dissipation rate, ε. A nonlinear Reynolds stress model (Shih, Zhu & Lumley 1996) is employed to relate the Reynolds stresses and the strain rates of the mean flow. To track free-surface movements, the volume of fluid (VOF) method is employed. To ensure the accuracy of each component of the numerical model, several steps have been taken to verify numerical solutions with either analytical solutions or experimental data. For non-breaking waves, very accurate results are obtained for a solitary wave propagating over a long distance in a constant depth. Good agreement between numerical results and experimental data has also been observed for shoaling and breaking cnoidal waves on a sloping beach in terms of free-surface profiles, mean velocities, and turbulent kinetic energy. Based on the numerical results, turbulence transport mechanisms under breaking waves are discussed.