In this paper, the hydrodynamic and free surface elevation fields in breaking waves are simulated by solving the integral and contravariant forms of the three-dimensional Navier–Stokes equations that are expressed in a generalized time-dependent curvilinear coordinate system, in which the vertical coordinate moves by following the free surface. A new k−l turbulence model in contravariant form is proposed; in this model, the mixing length, l, is defined as a function of the maximum water surface elevation variation. A new original numerical scheme is proposed. The main element of originality of the numerical scheme consists of the proposal of a new fifth-order reconstruction technique for the point values of the conserved variables on the cell face. This technique, named in the paper as WTENO, allows the choice procedure of the reconstruction polynomials for the point values to be modified in a dynamic way.