As part of an effort to develop a phase-resolving wave driver and establish a robust foundation for a comprehensive morphology model capable of describing the year-long circulation process of sandy beaches and addressing beach erosion, the authors introduced a wave driver comprising the spatially averaged Navier–Stokes equations. To verify the newly proposed wave driver, the authors numerically investigated the nonlinear shoaling characteristics of the surf and swash zones. The authors also tested the validity of the eddy viscosity model for Reynolds stress due to wave breaking using data from the Super Tank Laboratory Data Collection Project (Krauss et al., 1992). The characteristic length scale of the breaking-induced current is not negligible, posing a clear contradiction to the applicability of widely used eddy viscosity models, such as the, and instead favoring large eddy simulation (LES) with finer grids. In light of these observations, the residual stress in the spatially averaged Navier‒Stokes equation is modeled based on the Lagrangian dynamic Smagorinsky approach (Meneveau et al., 1996). The authors numerically integrate newly proposed wave driver using SPH with a Gaussian kernel function. A severely deformed free water surface profile, free-falling water particles from the wave crest, queuing splashes after water particles land on the free surface, and wave fingers resulting from structured vortices on the up-wave side of the wave crest (Narayanaswamy and Dalrymple, 2002) are successfully duplicated in the numerical simulation of wave propagation over a uniform slope beach: these features have been regarded as very difficult to duplicate in computational fluid mechanics. The numerical simulation also indicates that the widely used Standard Smagorinsky model with in the literature results in an excessively dampened water surface profile, attributed to the overestimated energy dissipation from wave breaking, leading to the loss of picturesque features, such as reverse breaking, observed both in nature and in numerical simulations using the Dynamic Smagorinsky model. Furthermore, the bottom shearing stress was estimated using the numerically simulated velocity profile and dynamic Smagorinsky coefficient, rather than relying on the quadratic friction law with a friction coefficient, as in the literature. The observation revealed that the maximum bottom shearing stress occurred when a broken wave, commonly known as a bore, rushed into the deep swash zone. Additionally, the study demonstrated that every aspect of the evolution of bottom shear stress within a wave period, such as its asymmetric nature over the surf zone where most of the sediment available along the beach is activated, could be precisely simulated using the newly proposed wave driver. These features of bottom shear stress over the surf and swash zones are crucial prerequisites for a morphology model to accurately describe the year-long circulation process of sandy beaches and effectively address beach erosion. This is particularly important, as the seasonal migration of an offshore bar is closely related to asymmetrically accelerated flows.