As wave propagates into shallow water, the shoaling effect leads to increase of wave height, and at a certain position, the wave will be breaking. The breaking wave is powerful agents for generating turbulence, which plays an important role in most of the fluid dynamical processes in the surf zone, so a proper numerical model for describing the turbulent effect is needed urgently. A numerical model is set up to simulate the wave breaking process, which consists of a free surface model using the surface marker method and the vertical two-dimensional model that solves the flow equations. The turbulence is described by Large Eddy Simulation (LES) method where the larger turbulent features are simulated by solving the flow equations, and the small-scale turbulence that is represented by a sub-grid model. A dynamic eddy viscosity sub-grid scale stress model has been used for the present simulation. The large eddy simulation model, which we presented in this paper, can be used to study the propagation of a solitary wave in constant water depth and the shoaling of a non-breaking solitary wave on a beach. To track free-surface movements, The TUMMAC method is employed. By applying the model to wave breaking problem in the surf zone, we found that these model results compared very well with experimental data. In addition, this model is able to reproduce the complicated flow phenomena, especially the plunging breaker.