A NEW 3D ROLLER APPROACH FOR FACING ROTATIONAL SURF ZONE HYDRODYNAMICS

A 2DH highly nonlinear Boussinesq-type of model for breaking waves has been developed in order to inve stigate surf zone hydrodynamics, also in the presence of complex bathymetries. The set of equations includes contin uity and rotational momentum equations, coupled with the vor ticity transport equation. An appropriate spatial d efinition of the 3D roller concept, along with an algorithm for accurately tracking the roller position, have been on purposely developed. Several numerical simulations have been carried out for the case of a submerged elliptic sh oal. The results have been compared with both experimental data and with the results of other numerical models availabl e in the literature. Finally, the vorticity dynamics under a breaking wave has been analyzed both in time and s pace, showing that a fairly correct interpretation of the wake ef fect in the rear part of the wave crest is obtained .


INTRODUCTION
Many coastal studies are aimed at investigating the complex interaction of the wave motion with non-uniform bed morphologies.Indeed, in such cases, phenomena such as refraction, diffraction, reflection, shoaling need to be taken into account.The situation is even more complex in the presence of breaking waves, since the flow becomes rotational and vortices are generated, thus inducing important dissipation of wave energy.The breaking process also induces a number of phenomena, such as wave set-up, rip current and alongshore current.In the presence of erodible bottom, a huge amount of sediment transport occurs as well.Just to provide a picture of the problem, a nonexhaustive overview of the aforementioned phenomena is sketched in Figure 1.A fairly good equilibrium between an accurate description of the above mentioned wave-forced phenomena in coastal regions and a limited computational effort is given by Boussinesq-type of models.Indeed they appropriately describe wave propagation and permit to extract information about the vertical distribution of the horizontal velocities despite they are managed as depth-integrated equations.After the first definition of such kind of equations by Boussinesq (1872), in the last decades their dispersive and nonlinear characteristics have been greatly improved, see for example: Madsen and Sørensen (1992), Wei et al. (1995); Gobbi and Kirby (1999); Lynett and Liu (2004); Madsen et al. (2006); Bingham et al. (2009).However such formulations do not include directly energy dissipation terms due to breaking phenomena because of the hypothesis of irrotational flow motion.Therefore, in order to face wave breaking effects, it is necessary to insert additional terms in the momentum equation able to simulate the enhancement of the momentum flux generated by turbulence due to breaking.Traditionally two main strategies have been adopted: 1. eddy viscosity models (Zelt, 1991;Kennedy et al., 2000) which introduce artificial viscosity terms into the momentum equation in order to model the turbulent mixing and the dissipation caused by breaking; 2. surface roller models (Schäffer et al., 1993;Madsen et al., 1997) which account for the excess of momentum originating from the non-uniform velocity distribution under a breaking wave.
The intent of obtaining a more physically based model led Veeramony and Svendsen (2000) to remove the hypothesis of irrotational motion and to solve the vorticity transport equation coupled with a weakly nonlinear Boussinesq one-dimensional model.In particular, such a model was derived starting from the Reynolds-averaged equations and the source of vorticity was assumed to be located on the turbulent front of the breaking wave.To this aim, the roller concept has been used as it was geometrically defined by Svendsen (1984), and the breaking generated horizontal vorticity has been injected through the lower edge of the roller itself inside the domain.The value of the vorticity at this location is quantitatively determined by assuming the hydraulic similarity between breaking waves and hydraulic jumps with low Froude numbers (Svendsen et al., 2000), while the vorticity transport equation is solved analytically under the assumption of a depth-constant eddy viscosity.Musumeci et al. (2005) extended such an approach by deriving a fully nonlinear version of the model, in order to obtain better shoaling properties close to the breaking point, and by implementing an original algorithm to avoid losses of vorticity at the toe of the roller.Although the results obtained by Veeramony and Svendsen (2000) and Musumeci et al. (2005) satisfactorily compare with laboratory data, real situations are quite different.Indeed in the field the angle of wave attack or the randomly variable water depth, due to natural or artificial complex bathymetries, dramatically affect wave propagation.Furthermore, breaking generated longshore currents are totally out of the picture of a one-dimensional model.Aiming at coping with such problems, in the present work a new set of 2DH horizontal equations has been derived, able to deal with the rotational motion which is characteristic of the surf zone.The equations include the effects of vorticity, by extending the 1D fully nonlinear approach of Musumeci et al. (2005).In such a framework, an appropriate numerical strategy has been also developed to track the position of the surface 3D roller over the horizontal plane.
In the following section, first a description of the model is provided, along with a discussion of the derived governing equations; then the numerical algorithm implemented to track the 3D roller and, in turn, to solve the vorticity transport equation is presented; the results obtained for the case of an elliptic submerged shoal are analyzed by making a comparison with both experimental and other numerical results; moreover the 3D vorticity dynamics under breaking waves is investigated; finally some of the main conclusions of the work are summarized.

Variable definition
The Cartesian coordinate system adopted is assumed in such a way that x ˆ and y ˆ are the horizontal coordinates, with x ˆ in the cross-shore direction pointing onshore, and z ˆ is the vertical one, being 0 ˆ= z at the still water level.Figure 2 shows the adopted reference coordinate system along with the water depth h ˆ and the main variables of the problem, i.e. the three velocity components u ˆ, v ˆ and w ˆ and the free surface elevation ζ ˆ.The mathematical notation adopted in the present work is such that the symbol ( ˆ ) means that the corresponding variable is dimensional.
The derivation of the proposed 2DH Boussinesq-type equations has been executed in dimensionless form.The following dimensionless parameters have been chosen to be representative of the investigated flow: (i) the dispersive parameter , where 0 ĥ is the water depth at a reference location, and 0 k and 0 â are the corresponding wave number and wave amplitude respectively.
In shallow waters it is usually assumed that the changes in horizontal dimensions, scaled by means of the wave number, are much slower than the changes in the vertical dimension, generally scaled by means of the water depth.This leads to the fact that it is reasonable to assume that µ 2 << 1.As regards the nonlinear parameter δ, whereas in the so-called weakly nonlinear model it is assumed as O(µ 2 ) = O(δ ), in the present model no restriction are introduced on the magnitude of δ.In this sense, the model can be classified as fully nonlinear.As demonstrated by Musumeci et al. (2005), the main advantage of such characteristics is related to the fact that wave propagation during shoaling close to the breaking point can be better predicted.The independent variables can be appropriately made dimensionless as follows: where g is the gravity acceleration and t ˆis time.The free surface elevation ζ ˆ, the pressure p ˆ and the three velocity components in dimensionless form read: where 0 â is the wave amplitude.Since in the present model the flow after breaking is considered to be rotational, the vorticity field ω ˆ must also be solved.From the definition of vorticity, its components x ω ˆ, y ω ˆ and z ω ˆ can be scaled as follows In order to model the turbulent shear stresses, an eddy viscosity t νˆ is introduced, which in dimensionless term is: where C ν is a dimensionless coefficient which can be estimated to be about 0.01 ÷ 0.03, according to the experimental analyses of Cox et al. (1995).
In order to derive the Boussinesq equations, a law for the velocity profile has to be considered, depending on some characteristic velocity.Several choices are possible, as the depth-averaged velocity or the velocity at some specific elevation.In the present case the depth-averaged velocity, whose components along x and y are: has been chosen.
The corresponding depth-averaged horizontal velocity vector is then . Moreover, by considering that the flow is rotational, the dependency on the vorticity is made explicit through the horizontal rotational velocity profile vector ) . Such an horizontal rotational velocity turns out to be a component of the total horizontal velocity.In particular, eq. ( 9) has been derived here in a more general and flexible form compared to the corresponding procedure proposed by Musumeci et al. (2005).Indeed, in that work the authors adopted the stream function approach which is strictly two dimensional on the vertical plane and it cannot be extended to the present 2DH case.The method followed here consists: (i) in expressing the vertical component of velocity in terms of the horizontal ones by applying the continuity equation; (ii) in submitting the obtained expression of the vertical velocity into the definition of the horizontal vorticity components, i.e. eqs.( 4)-( 5); (iii) in integrating over the depth between the bottom -h and the generic elevation z.

The formulation of the problem
To derive the proposed rotational Boussinesq-type of formulation, the continuity and the Reynolds averaged Navier Stokes equations are integrated over the depth.Therefore, by applying the assumptions of fixed bed with a free slip condition and the Leibniz rule at the surface and at the bottom, the pressure term is eliminated from the equations.Moreover the expression for the horizontal velocity profile is used in order to obtain the governing equations as a function of the depth averaged velocity, of the vorticity and of the free surface elevation.After some algebra, the continuity and the combined momentum equations read: where the pedix ( t ) indicates time derivatives.Moreover, in order to improve the dispersion characteristics of the model in deeper water, the linear operator Madsen and Sørensen (1992), has been adopted.
It must be noticed that in eq. ( 11) the dependency on the rotational velocity r u , which in turn is function of the vorticity ω r injected inside the flow by the breaking mechanism, is included exclusively within the terms and uw D , which, in analogy to the work of Veeramony and Svendsen (2000), are called breaking terms.Indeed, they represent the excess of momentum flux due to breaking, which in turn is related to the dissipation of wave energy in the surf zone.In particular: give the excess of momentum flux due to the vertical variation of the rotational velocity; is the contribution to the pressure due to the vertical motion; w D is the excess of momentum due to the vertical motion; sv D and sh D are the shear stresses inside the fluid which depend, respectively, on the turbulence along the vertical and the horizontal planes; uw D represents the interaction between the waves and the mean flow.
Considering the solution of the vorticity transport equation, by assuming that within the surf zone the effect of convection is leading order with respect to the effect of vorticity stretching, the problem of determining the two-dimensional vorticity field can be reduced to a one-dimensional problem along the curvilinear abscissa s, which corresponds to the direction of wave propagation throughout the domain (see Figure 3).By cross-differenciating the Reynolds-averaged Navier-Stokes equations and by applying the discussed scaling arguments, the dimensionless vorticity transport equation along the abscissa s reads: where u s is the dimensionless horizontal velocity in the direction of wave propagation, which has a module equal to that of the horizontal velocity u.
An analytical solution of such a vorticity transport equation has been found through a perturbation approach.The solution up to O(δ ) has been obtained starting from the basic state, which corresponds to

Direction of wave propagation
the linear solution.In this way, the vorticity distribution is calculated at a fixed time step, starting from the vorticity field at the previous steps and from the vorticity ω s pumped into the water column from the roller region by the breaking phenomenon.In order to find the solution, the vertical coordinate is stretched, in such a way to follow the shape of the domain, i.e. of the breaking wave.Finally the solution is found in a Fourier series form.

THE NUMERICAL SOLUTION
The numerical scheme chosen to integrate the momentum and the continuity equations is always a crucial point; indeed it can make the difference between an accurate and efficient model and a poor one.In Boussinesq-type of models, a widely adopted numerical solution strategy (see, for example, Wei et al., 1995) is the one which considers a fifth-order finite difference scheme for the spatial derivatives and a thrid-order predictor, fourth-order corrector scheme, for the integration in time, known in the literature as Adams-Bashforth-Moulton or ABM scheme.The same approach has been used here since it allows to take advantage both of the relatively fast explicit scheme and of the accuracy of the spatial derivative discretization.
As regards the offshore and the lateral boundary conditions, the computational grid is a rectangle in which three of the boundaries, namely the lateral ones and the onshore one, are vertical walls and the forth one is the offshore boundary.
The condition at the vertical walls is that the horizontal velocity component is zero, which is an artificial condition for many applications.In order to overcome such a problem a sponge layer can be inserted in front of each wall.At the offshore boundary an absorbing-generating boundary condition is applied, modified with respect the one proposed by van Dongeren and Svendsen (1997) for their nonlinear shallow water equation model, since, in the present case, the forcing term includes also dispersive effects.

The 3D Roller
As mentioned, the wave breaking phenomenon has been modelled here by trying to follow its actual physical behaviour.In particular the loss of wave energy has been assumed to be related to the vorticity which appears below the free surface during the breaking of the waves.The main assumption of this model is that the vorticity is "pumped" inside the water column from the highly turbulent region which is present on the top of the wave during the breaking phenomenon.
Such an approach has been already followed by Veeramony and Svendsen (1999) and by Musumeci et al. (2005) for the one horizontal dimensional case.In those models the direction of wave propagation is known and the turbulence related to breaking has been modeled by a single twodimensional roller, in which the vorticity ω s is perpendicular to the vertical plane of the roller.
In the proposed model the presence of both the horizontal directions (x, y) causes a modification of the roller itself, which becomes a mass on the top of the wave which can change its form and length in time, since the same wave can start to break at different time and space along its front.The approach followed here to overcome this problem is to consider the so called "3D roller ", as composed of several two-dimensional rollers which are placed on vertical surfaces which follows the wave propagation direction (see Figure 4).
Moreover, the implementation of the analitical solution of the vorticity transport equation (i.e.eq.12) for a generic two horizontal dimensions case needs an appropriate spatial definition of the 3D roller in order to determine the value of ω s to be specified along the upper limit of the computational domain, corresponding to the lower boundary of the surface roller itself.
A critical point, which is common to all the Boussinesq-type of models, is to establish a criterion for the onset of breaking, which is here complicated by the need to determine also the spatial location of the 3D roller, and particularly of the toe of the roller.In the proposed model it is assumed to trigger breaking when and where the local steepness of the wave front exceeds a critical surface gradient, equal to φ tan .Precisely, the toe of the roller is defined as the location where the wave steepness is identical to φ tan and the roller itself is included between its toe and the wave crest.In order to account for the transition from the initial breaker to the bore-like stage within the inner surf zone, the critical roller angle, φ , is assumed to decrease from B φ up to its final value 0 φ as a function of the age τ of the roller, i.e. the time interval measured from the initiation of wave breaking.The relation used here, initially introduced by Schäffer et al. (1993), is ( ) where t 1/2 defines the time scale for the development of the roller, and t B is the time of incipient breaking; the values for the coefficients presented in eq. ( 13), which have been used in all the results presented here, are: φ B = 30° , φ 0 = 12° and t 1/2 = T/5, where T is a typical period of the waves.
It is worth pointing out that such a criterion permits a reliable definition of the toe position, i.e the zone where the injection of vorticity is stronger.Therefore, notwithstanding its complex numerical implementation over two horizontal dimensions, in the present case such a breaking criterion has been preferred to other criteria, as the one based on the time derivative of free surface elevation (see Kennedy et al. 2000).
Indeed, in the two horizontal dimensions, the toe of the roller becomes a curve (see Figure 4), which can be defined as the points satisfying the condition that the absolute value of the surface elevation gradient | ζ ∇ | is equal the instantaneous local value of φ tan and that the gradient in the direction of the wave propagation is negative.This last assumption is related to the fact that the present model aims at simulating the physical behaviour of spilling breakers, for which the roller appears on the front of the wave, where the free surface elevation decrease in the direction of propagation, i.e if 0 / < ∂ ∂ s Another characteristic of the roller toe angle φ is that it can vary within the same roller during the different stages of the breaking process.For example, when oblique regular waves break on a beach with straight and parallel depth contours, both the initial and the final stages of the breaking process may be represented within the same roller.A clear example of such a situation is shown in Figure 4. Hence, it is necessary to keep track of the surface roller age along the orthogonal to the roller toe curve.Moreover the evolution of the roller position depends on the time-varying wave shape; thus an algorithm for accurately tracking the roller characteristics has been developed here.Even though the surface roller detection appears quite simple to be determined by using a continuous description, such a detection becomes rather complex over a numerical grid, where a discrete representation of the spatial variation of the hydrodynamic variables along the roller is required.
In the past Sørensen et al. (1998) faced a similar problem by using a staggered grid for the determination of the toe points, by assuming a linear variation between two toe points and by defining the roller as the area comprised between each couple of toe points and the projection on the horizontal plane of the two lines tangential to the direction of wave propagation.
In the proposed model, starting from the above mentioned strategy, the following approach has been adopted.With reference to Figure 5, it must be preliminarly stressed that the Boussinesq equations are derived over a rectangular grid; therefore the dependent variables ζ and u are known at each grid point.Instead, the surface elevation gradient, i.The algorithm for accurately tracking the roller position has been developed and implemented as follows: 1. a rectangular grid is defined, with a surface elevation computational node at each grid node; 4. the wave crest and the width of the roller in the wave propagation direction α are also found; 5. the injection of vorticity ω s is computed along the sections comprised between the toe points and the wave crest; 6. the age of the roller τ is updated and the corresponding toe angle can be calculated at the next time step.
Finally, it must be stressed that the variable ω s is computed along the roller subgrid, but it is transferred and stored in the correspondence of the rectangular fixed grid, in order to solve the vorticity transport equation and therefore in order to obtain the dissipative breaking terms of the proposed Boussinesq-type of equations.

WAVE PROPAGATION OVER A SUBMERGED SHOAL
To test the capabilities of the proposed model to handle a fully 2D bathymetry, numerical simulations have been carried out by considering a submerged elliptic shoal, as in the experiments of Vincent and Briggs (1989), both in breaking and non breaking conditions.Indeed, this test is an ideal case for the validation because intrinsic two-dimensional phenomena such as diffraction, refraction and jet-like strong shear current induced by breaking, are generated.Moreover other authors have also adopted such a test as a benchmark for validating numerical models.As an example, Choi et al. (2009) performed numerical simulations using both phase-averaged literature models and phase-resolving literature models.In particular they used also the Boussinesq-type of model developed by Chen et al. (2000), known in the coastal community as FUNWAVE.
Figure 6 shows the top view of the bathymetry used here for the application of the proposed model.Moreover in the mentioned Figure 6 a section onshore of the shoal is shown (Transect 4), at which experimental data are provided.For the numerical simulation, the adopted computational grid is uniform (0.1 m x 0.1 m) and the domain where the results have been analyzed is 19 m along the x ˆ direction and 25 m in the y ˆ direction.As a matter of fact, in order to avoid distortions due to the boundary the actual numerical domain is 14 m larger in the y ˆ direction and 20 m larger in the x direction.Moreover the time step for the numerical integration is fixed equal to 0.013s .
Two different regular waves have been considered as input of the proposed model: the first one, which represents a non breaking case, has significant wave height equals to 2.54 cm, wave period of 1.3 s and corresponds to case M2 of the experiments; the second wave, which corresponds to the case M3 of the experiments, has the same wave period as the previous one but the wave height is much larger and equal to 13.5 cm thus providing a test for the breaking model.

Wave height distribution
In order to give a comprehensive description of the wave height modification due to the presence of the submerged elliptic shoal, the results of the proposed numerical model have been provided in terms of the disturb coefficient C d , i.e. the ratio of the local depth over the incident wave height.Figure 7 shows the top view of such an outcome, for the non breaking M2 case, obtained by using the proposed model and the results obtained by Choi et. al (2009) by using the FUNWAVE model (Chen et al., 2000).
It can be noticed that the presence of the submerged shoal causes reflection-refraction effects which are highly dominated by the shape of the shoal.Indeed the regions which have a disturb coefficient greater than 1 are strongly bent around the shoal.Shoreward, the presence of the shoal induces a refraction process which generates a strong symmetric focusing of wave energy after the crest of the shoal, whose values are in agreement with the results of Vincent and Briggs (1989).Moreover, as it can

Direction of wave propagation
Absorbing-generating boundary condition be seen, the interference patterns are stretched and oriented towards the centre of the shoal.Such results are very similar to the numerical results obtained by Choi et al. (2009) by using the FUNWAVE model.
In Figure 8 numerical results of the proposed model in terms of the disturb coefficient C d have been compared also with the experimental results gathered along Transect 4 (by Vincent and Briggs, 1989).Both numerical and experimental data show the presence of an increase of the wave height shoreward of the submerged shoal.Generally, small values of the gradient of the disturb coefficient are obtained by the proposed model compared to the laboratory data.Indeed, for example, the location of the lowest values of the wave height dislocated of about towards the lateral boundaries of the domain.
The analysis of the breaking wave case M3 is provided in Figure 9 where the results of both the proposed model (a) and the FUNWAVE model (b) are showed.
In such a case, the breaking effect causes a completely different pattern shoreward of the shoal in comparison to the non-breaking one.First of all, it can be observed that offshore of the shoal, while the results of Choi et al (2009) in Figure 9.b show that the wave fronts are not bending, thus obtaining that reflection does not act when waves are breaking, in the present case (see Figure 9 preserved.As a matter of fact it is questionable that breaking may eliminate entirely sea-ward effects, as shown by the results of Choi et al. (2009).An explanation of such a result could be attributed to the eddy-viscosity breaking-generated model used in FUNWAVE which spreads the breaking dissipation more than it should do through the entire domain.Such a behavior is not present in the proposed model, since the introduction of breaking generated vorticity is strictly limited backward from the wave crest position.
Furthermore, in such a case, the shoaling effect causes the wave breaking condition to be satisfied just before the top of the shoal.After that, the breaking phenomenon evolves and the excess of momentum flux due to it causes a jet-like current on-shore directed.Such a current contrasts the convergence of the wave rays within the focusing zone, seen in the non breaking case.As a result, a decrease of wave height at the back of the shoal is obtained.Such a mechanism was explained by Yoon et al. (2004) and it is caught also by the present model, as it can be observed by the results in Figure 9.a.
It can be also noticed that the breaking waves are clearly strongly three-dimensional because of reflection and refraction induced by the depth varying bathymetry.It follows that the numerical model, and in particular the roller tracking module, seems to be able to represent such a case of a curved wave front.Moreover the horizontal distribution of the normalized wave height shows that an increase of wave height is located more shoreward in comparison to the non breaking case, along with the alternation of region with lower and higher waves.
The comparison of the proposed model results with the available experimental data along the Transect 4 for such a breaking case M3 is shown in Figure 10.From its analisys it is possible to state that the proposed model can interpret the order of magnitude of the wave height after breaking, although it appears that the dissipation induced by wave breaking is still not high enough to reproduce the experimental data accurately.

Vorticity results
The breaking case M3 of Vincent and Briggs (1989) represents a good reference for testing the capabilities of the proposed model to simulate realistic vorticity dynamics induced by the injection of turbulence inside the flow domain.Indeed, in such a case the direction of wave propagation is deviated by refraction caused by the shoal, which also determines an increase of wave steepness up to breaking.In order to carry out a detailed analysis, the attention has been focused on a water column located onshore of the crest of the shoal (see Figure 11).The contour plots of the spatial distribution of vorticity ω at the still water level (i.e. at z = 0) within such a region are shown in Figure 12, along with the corresponding contour plots of the surface elevations.The attention has been fucused over half of the wave period (T ), in particular it has been considered the time interval during which breaking generated vorticity is espected to be not negligible, i.e. from the instant of incipient breaking (t = 0.2 T ).Indeed, the vorticity field at such a timestep (Figure 12.f) is characterized by low values, due to the fact that vorticity at the roller have not yet spread into the water column.During the following instants (from 0.3 T up to 0.5 T ) the waves are still breaking, as demonstrated by the decreased elevation of the wave crest, while the vorticity increases suddenly, reaching its maximum (40 Hz).Such a behaviour shows that the vorticity dynamics is solved appropriately indeed, in actual breaking waves, the presence of vortices on the rear part of the wave appears a bit delayied with respect to the initiation of the breaking phenomenon and the foam generated vorticity increases gradually also at the front of the wave.A demostration of this peculiarity can be observed in Figure 4, where different stages of breaking phenomenon can be observed in the same 3D wave.During the second half of the wave period, i.e. for t greater or equal to 0.5 T , the slope of the free surface elevation is not large enough for maintaining the breaker and therefore wave breaking stops.As a consequence, vorticity decreases both in intensity and in space.
As regards the analysis of vorticity spatial evolution, it is possible to note that vorticity is present mainly in the correspondence of the wave front at the initial instants of breaking (see for example Figure 12.g for t = 0.3 T ).During the following time steps, up to t = 0.4 T, the planar extension of vorticity increases in both the horizontal directions, penetrating also offshore with respect to the wave crest.Such a behaviour seems to support the fact that the proposed model is able to simulates the wake generated on the rear side of surf zone waves.

CONCLUSIONS
A new two-dimensional Boussinesq type of model has been derived and implemented, aimed at studying the flow driven by breaking and non breaking waves in the surf zone, also in presence of complex bathymetries.The main feature of such a model is the absence of the usually adopted limiting hypothesis of irrotational flow; thus a more physical description of the flow within the nearshore region is provided.To this aim, particular attention has been devoted to the description of the breaking process, here implemented by using the well known surface roller approach.In particular, an extended 3D roller model has been developed and implemented by means of an original numerical strategy in order to accurately compute the vorticity injected inside the domain by the breaking phenomenon.
The importance of applying the above mentioned approach to a two-dimensional Boussinesq-type of model, as in the case being, is that it can manage in a proper way the presence of macrovortices generated by breaking waves over highly three-dimensional bathymetries like reefs, submerged breakwaters, gulfs, cuspidate beaches, etc.
The proposed model has been tested for waves propagating over a submerged shoal for both non breaking and breaking conditions.Precisely the experiments of Vincent and Briggs (1989) and the numerical results obtained by Choi at al. (2009) using the well-established FUNWAVE model have been considered for comparison.From the analysis of the results on the spatial distribution of wave height, expressed in dimensionless form as the disturb coefficient C d , it is possible to state that the numerical data are in a fairly good agreement with the experiments, above all for what concerns the presence of oscillations between higher and lower wave height regions shoreward of the submerged shoal.
The vorticity dynamics under breaking waves has been deeply investigated for the wave breaking shoal case M3.In particular, the attention has been focused over a water column located onshore with respect to the top of the shoal where wave breaking develops.The analysis of the spatial and temporal evolution of the computed vorticity has shown that the vorticity increases rapidly in a small region near the roller and immediately below the free surface.Moreover it has been found that the horizontal extension of the vorticity structures increases, thus simulating the wake effect in the rear part of the wave crest.
Notwithstanding some improvements which could be introduced in the numerical solution of breaking terms, the application of the proposed model to study nearshore hydrodynamics generated by complex conditions of engineering interest, even in the presence of breaking, seems to be quite satisfactory.

Figure 1 .
Figure 1.Sketch of the main nearshore dynamics phenomena induced by waves.

Figure 2 .
Figure 2. Definition of the main variables and reference system adopted in the present model.

Figure 3 .
Figure 3. Plan view of a wave ray field with indication of the curvilinear abscissa s and of the horizontal velocity us.

Figure 4 .
Figure 4. Photo of a breaking three-dimensional wave, with a qualitative indication of the 3D roller, i.e. the region placed between the wave crest and the roller toe curve; c is the wave celerity, φ φ φ φ is the angle at the toe of the roller which influences the breaking onset criterion; ω ω ω ωs is the vorticity injected at the roller lower boundary.
the correspondence of the central nodes of the rectangular grid (i.e. at the circles in Figure5).

Figure 5 .
Figure 5. Top view of the numerical grid in which the roller toe points (×), rectangular grid points (•) and centers of grid cells (o) are showed; the broken lines represent the sections at which the injection of vorticity ω ω ω ωs is computed.
2. values of the roller angle φ and of the surface elevation gradient ζ ∇ are then calculated at the centre of the grid cells; 3. roller toe points satisfying the toe condition φ

Figure 6 .
Figure 6.Bathymetry of the submerged elliptic shoal defined by Vincent and Briggs (1989) with indication of the Transect 4 at which results are provided.

Figure 7 .
Figure 7. Non-dimensional wave heights (i.e.disturb coefficient Cd) for the non breaking wave case M2, with indication of the section at which the experimental results of Vincent and Briggs (1989) are provided: (a) proposed model; (b) FUNWAVE model applied by Choi et al. (2009).

Figure 8 .
Figure 8. Wave height distribution normalized respect to the incident one (Cd) at Transect 4 for the M2 wave.
.a) the curvature is

Figure 9 .
Figure 9. Non-dimensional wave heights (i.e.disturb coefficient Cd) for the breaking wave case M3, with indication of the section at which the experimental results of Vincent and Briggs (1989) are provided: (a) proposed model; (b) FUNWAVE model applied by Choi et al. (2009).

Figure 10 .
Figure 10.Wave height distribution over the incident one (Cd) at Transect 4 for the M3 breaking wave.

Figure 11 .
Figure 11.Top view of the shoal with indication of the water column in which the description of vorticity is investigated.

Figure 12 .
Figure 12.Distribution of the free surface elevation ζ ζ ζ ζ (form a to e) and of the vorticity ω ω ω ω (from f to j) within half a wave period, i.e. from 0.2 T to 0.6 T, as obtained at the still water level (z = 0) of the water column.