The breaking of gravity water waves induces both strong turbulence near the free surface and air–water mixing, which are not captured in currently-available single-phase, non-hydrostatic Reynolds Averaged Navier Stokes (RANS) models. In order to account for such dynamics, the boundary conditions proposed by Brocchini (2002) have been implemented in a two-dimensional vertical (2DV) non-hydrostatic RANS numerical model. For the theoretical boundary conditions of Brocchini (2002) to be applicable in the numerical model, appropriate closures have been adopted and mathematical equations have been proposed to calculate the coefficients introduced in the mentioned boundary conditions. Navier Stokes equations along with different turbulence closure models have been solved using finite volume method and pressure correction technique. According to the new boundary conditions, the normal-to-mean surface gradient of the turbulent kinetic energy (TKE) differs from zero and is computed as a balance between production and dissipation of TKE within the air–water mixing layer. Also, the kinematic and dynamic boundary conditions have been modified accordingly, to account for the effects of the thin two-phase layer formed at the free surface. The modified kinematic boundary condition allows for the mass exchange between the two-phase layer and the main body of the water. For the first time it is demonstrated that the proposed simple analytical model of Brocchini (2002) leads to improvements in the prediction of incipient breaking and the wave characteristics in the surf zone. This opens the way to a new, accurate but numerically low cost, approach for the computation of the air–water mixing that characterizes breaking waves. All the many benchmarking tests run to verify the ability of the new model show that significant improvements are achieved. Improvements are specifically observed in the prediction of: (I) breaking point and the breaking wave height, (II) the consequent dissipation of the wave energy observed in the form of the crest level distribution throughout the surf zone, (III) the magnitude of the horizontal velocity near the free surface, and (IV) the TKE distribution in depth. Notably, the novel numerical model does not use any parametric criteria for detecting the inception of breaking, thus it does not need calibration for different flow circumstances. Also, since no special treatments, such as hydrostatic pressure assumptions at the front face of the breaking waves are considered, the new model well captures the breaking-induced dissipation while giving a more accurate estimation of the dynamic pressure.