A Comparative Study of Eddy Viscosity Type Models for Simulating Wave Breaking Induced Energy Dissipation in Boussinesq Models

Boussinesq-type models in their original form cannot be used to simulate wave breaking or wave breaking induced energy dissipation across the surf zone, as they are developed under the hypothesis of ideal and irrotational flow. Since these equations restrict the production and the evolution of vorticity within, researchers have used a few closure models to introduce depth limited wave breaking, which primarily incorporates friction type, surface roller type or eddy viscosity type. Within those, a couple of eddy viscosity type models have been used successfully to introduce wave breaking induced energy dissipation nearshore over the years, but with less emphasis on the effect of bottom configuration on wave evolution. Hence, the applicability of two eddy viscosity type closure models for simulating wave breaking induced energy dissipation over seabeds with plane slopes and submerged structures is investigated coupling them into two widely used Boussinesq-type equations, in the present study. A published and a new set of data are used to perform the comparative study of these two sub-models for one dimensional horizontal (1DH) wave transformation.


Introduction
Accurate prediction of wave climate in a nearshore region is of prime importance to coastal engineers as the evolution of nearshore waves significantly controls the dynamic equilibrium of a given site. Long term wave data are mostly unavailable for nearshore sites. Therefore, researchers and engineers generally resort to deepwater wave data to predict nearshore waves and currents. In order to perform such predictions, numerical models are used extensively but with limitations. A good numerical model for a nearshore region must be able to incorporate effects of wind-wave generation, nonlinear wave-wave interactions, shoaling, refraction, reflection, diffraction, wave breaking, energy dissipations due to wave breaking and bottom friction, wave-current interactions and wave runup. However, the existing wave/wave-current models (e.g. based on mild slope equations, energy balance equation, Boussinesq equations, RANS equations, etc.) are developed under different assumptions, and none of them can accommodate all transformation processes in their original form.
Starting with the classical form of Boussinesqtype equations derived by Peregrine [1], researchers have made extensive developments and improvements to these equations over the last few decades. The poor dispersion characteristics were the major limitations of the original form of the Boussinesq equations for wave modelling, but Madsen et al. [2] and Nwogu [3] with two different approaches extended the validity of these equations from shallow water to deep water. Most of the subsequent theories, which were developed to extend the nonlinearity and dispersion characteristics { [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14]}, were based on Madsen et al. [2] and Nwogu [3] platform equations.
Progressing beyond the modelling of wave evolution over impermeable profiles, incorporation of porous flow resistance in Boussinesq-type models has been attempted and successfully introduced by Cruz et al. [15], Hsiao et al. [16], Garcia et al. [17], Chen [18], and Cruz and Chen [19]. This has enabled the application of such models for simulating wave propagation over permeable beds and porous submerged structures even.
The models developed based on Boussinesqtype equations are capable of simulating all the wave/wave-current processes except wave breaking, and energy dissipation due to wave breaking and bottom friction in its original form. However, to overcome such limitations, researchers have successfully introduced a few closure models. Regardless of the formulation, each of the approaches can be thought of as a means of adding the breaking wave force (momentum mixing) terms to the momentum equations. All these models differ on how they treat the onset and cessation of wave breaking and the rate of energy dissipation [20]. With the growth of these closure models, many researchers have developed a number of nearshore wave-current models for research and commercial purposes (FUNWAVE, BOUSS-2D, MIKE21-BW, COULWAVE, etc.). However, very few studies can be found on a comparative study of the performances of these models for different bottom configurations. An exception is the study by Siddique et al. [21] who investigated the difference in wave transformation and hydrodynamic characteristics of three different wave breaking sub-models (i.e. friction type, eddy viscosity type and surface roller type) coupled with the weakly nonlinear Boussinesq equations of Madsen and Sorensen [10].
Zelt [22] first proposed a sub-model for wave breaking and wave breaking induced energy dissipation, in his study of runup of nonbreaking and breaking solitary waves on plane impermeable beaches with a Lagrangian finiteelement Boussinesq wave model. Following Zelt [22], Kennedy et al. [23] proposed a new eddy viscosity formulation but with extensions to provide a more realistic description of the onset and cessation of breaking. A strong analogy can be found between Kennedy et al. [23] model and Shaffer et al. [24] surface roller model. Following Boussinesq eddy viscosity concept, Nwogu [25] approximated the turbulence shear stresses to be proportional to velocity gradients in vertical direction (proportionality coefficient being the eddy viscosity) and introduced a term representing a vertical gradient in shear stresses due to breaking induced velocity fluctuations in the momentum conservation equations. In his study, Nwogu [25] related eddy viscosity to turbulent kinetic energy and a length scale (one equation turbulence closure model).
In the present study, the wave evolution is simulated using a truncated version of Chen [18] or Nwogu [3] Boussinesq-type wave equations coupled with Nwogu [25] type eddy viscosity model (one-equation turbulence model) or Kennedy et al. [23] type eddy viscosity model extending them to nearshore. The predictive skills of these wave breaking models are then compared using a published and a new set of data in one-dimensional horizontal (1DH) wave propagation. The data include wave height, mean water level distributions and time series of free surface variations over seabeds with plane slopes and impervious/pervious submerged structures.

2.1
The Governing Equations A truncated version of Chen [18] and Nwogu [3] Boussinesq-type equations are used as platform equations coupled with Nwogu [25] or Kennedy et al. [23] eddy viscosity-type sub-model for wave breaking in the present study. The truncated Chen [18] equations are analogous to Nwogu [3] equations when porous damping terms are switched off. A short description of the governing equations used in the two main models (Boussinesq-type equations) is presented below.
In one horizontal dimension (1DH), the truncated version of Chen [18] equations (i.e. by dropping higher-order nonlinear and porous damping terms) read as: where, is the porosity of the porous material, is the free water surface elevation, ℎ is the clear water depth, ℎ is the porous layer depth, ℎ = (ℎ + ℎ ) is the total water depth, is the water particle velocity at reference water depth , and is the water particle velocity at reference water depth inside porous layer along the wave propagation direction. The factors and are determined by matching dispersion characteristics of wave motion in both clear water and porous layer (not discussed in this paper).
is the porous resistance force evaluated at = as follows; where and are respectively, the linear and nonlinear porous resistance coefficients, and is the added mass coefficient. The porous resistance coefficients are estimated from the following relationships following Sollitt & Cross [26]: where is the kinematic viscosity of water (∽ 10 −6 m 2 /s). The intrinsic permeability, , added mass coefficient, , and nonlinear drag coefficient, , are determined by empirical formulae proposed by van Gent [27].
where, 50 is the median diameter of the porous material, , and are empirical coefficients that, in principle, are determined by experiments. However, van Gent [27] concluded that, although the coefficients and depend on the parameters like shape, aspect ratio or the orientation of the porous material, the values 1000 and 1.1 can be used for and respectively, if the characteristic length scale 50 is used to define the size of porous material. In addition, he suggested using a constant value for empirical coefficient, of 0.34 for simplicity neglecting the complex dependency of coefficient, on the flow field. The porosity, n, of the porous material used in the present study, was found to be 0.44.
The above equations yield to Nwogu [3] Boussinesq-type equations when the porous damping terms are ignored. i.e.,

Energy Dissipation
Two additional terms are introduced into momentum conservation equation {Eq. (10) when Chen [18] type is used and Eq. (20) when Nwogu [3] type is used} to simulate the dissipation of wave energy due to bottom friction and wave breaking, which are represented by and . where, is the wave friction factor, and is the near-bottom wave orbital velocity along the xaxis. Following Kennedy et al. [23], a simple eddy viscosity type of formulation is adopted to simulate energy dissipation due to wave breaking (by introducing a momentum-mixing term into momentum conservation equation in the direction).

Eddy viscosity sub-model of Nwogu [25]
In Nwogu [25] type energy dissipation submodel, the magnitude of the eddy viscosity, , is related to the turbulent kinetic energy produced by wave breaking, , and a turbulent length scale, . The turbulent kinetic energy is determined from a semi-empirical transport equation with a source term for turbulent kinetic energy production by wave breaking, i.e., one equation turbulence closure model or k-model. This sub-model will be called TKE model hereafter.
where, is the free surface wave orbital velocity along the x-axis.
The parameter is introduced to ensure that the turbulence is produced only when exceeds a threshold velocity (a fraction of the wave celerity, C). The empirical coefficients, , and for this sub-model are set to be 0.08 and 1, respectively. To trigger the turbulent kinetic energy equation, an initial estimate of eddy viscosity, , is necessary. This is obtained by assuming a local balance between production and dissipation of turbulent kinetic energy as follows: .. (25) The wave breaking detection parameter, / , and dimensionless turbulent length scale, / are set to be equal to 1.0 and 1.0, respectively, irrespective of the configuration of the sea bottom following BOUSS-2D model {Nwogu and Demirbilek [28]} in all simulations.

Eddy viscosity sub-model of Kennedy et al. [23]
In Kennedy et al. [23] eddy viscosity sub-model, is determined in a similar manner to Zelt [22] but with some differences.
Where is a mixing length coefficient, which is found to be equal to 1.2. The quantity B varies smoothly from 0 to 1 to avoid an impulsive start of breaking and possible instability. Therefore B is defined by, The parameter * determines the onset and cessation of wave breaking. In other words, a breaking event starts when exceeds some threshold value, but as breaking develops, the waves will continue to break even if drops below this value. Therefore, the magnitude of * decreases in time from initial value to a terminal quantity . This gives, .. (28) where * is the transition time, 0 is the time that breaking was initiated and − 0 is the age of the breaking event. Kennedy et al. [23] proposed, = 0.65√ ℎ, = 0.15√ ℎ and * = 5√ℎ/ using their preliminary experimental data obtained over plane slopes. Chen et al. [29] recommended = 0.65√ ℎ for plane seabed slopes and = 0.35√ ℎ for barred profiles incorporating the effects of different breaker type. Choi et al. [30] re-calibrated Kennedy et al. [23] values and suggested that setting = 0.45√ ℎ would make the predictions better. Following Kennedy et al. [23], Chen et al. [29], Choi et al. [30] and by running the model multiple times with the recommended values for respective parameters, is set at 0.45√ ℎ, and 0.35√ ℎ for seabeds with plane slopes and with submerged structures, respectively, in the present study. and * are set at 0.15√ ℎ and 5√ℎ/ , respectively for any type of seabed configuration. Kennedy et al. [23] eddy viscosity sub-model is called Kennedy model for convenience hereafter.

Offshore Open Boundary and Moving Shoreline Boundary
As there is significant wave reflection from submerged structures, the waves which propagate back to wave incident boundary are absorbed to prevent a build-up of wave energy inside the computational domain using line boundary method proposed by Ishii et al. [31]. Kennedy et al. [23] slot technique is employed to simulate moving shoreline.

2.4
Numerical Scheme The truncated Chen [18] Boussinesq-type equations are solved in the time domain using a third-order Adams-Bashforth predictor step and a fourth-order Adams-Moulton corrector step. The computational domain is discretized with grid size, Δ in the -direction. The equation variables , ℎ, and are defined at grid points in a staggered scheme. The water depth and surface elevation are defined at the same grid points, while velocities are defined half a grid point on either side of elevation grid points. Therefore, the external boundaries of the computational domain correspond to velocity grid points. The first-order spatial derivatives are differenced to (Δ 4 ), which automatically eliminates error terms that would be of the same form as dispersive terms (Wei & Kirby, 1995). The second-order spatial derivatives are discretized to (Δ 4 ), while the advection terms are differenced with second-order upwind difference scheme.
The computations are carried out for free surface elevation, , reference velocity for clear water, , and reference velocity for porous medium, (cases with porous structures only). In addition, computations are carried out in a similar manner for turbulent kinetic energy, k, simultaneously when the TKE model is incorporated.

Experimental Set-up
The predictive skills of the models are compared using a published and a new set of experimental data, which include 07, 03, and 05 cases with plane slopes, solid submerged structures and porous submerged structures, respectively. Solid submerged structure models were constructed with concrete, and porous ones were constructed with gravel of mean diameter (d50) equal to 12 mm. Figure 1 illustrates the definition of variables associated with the physical models used and Table 1 lists the wave parameters, physical model dimensions and numerical space and time steps associated with all cases considered. All tests were conducted with regular waves.    Table 2.
The Willmott index for wave heights is found to be less than 0.95 only in case PLS07 when simulations are performed with TKE model, whereas this index is less than 0.95 in PLS04, PLS05, PLS06 and PLS07 when simulations are performed with Kennedy model. However, the performance of the Kennedy model in simulating mean water level is better than the TKE model for seabeds with plane slopes, in general. Both models, however, fail to predict the mean water level distribution across the surf zone for case PLS 06, indicating the necessity to improve the predictive skills of the models by fine-tuning the model parameters. Overall performance of the Kennedy model is found to   Table 2.
The Willmott index for wave heights is found to be slightly less than 0.95 in cases SSB02, SSB03 and PSB01 when simulations are performed with TKE model, whereas it is less than 0.95 only in SSB03 with Kennedy model. However, the performance of the Kennedy model in simulating mean water level is poorer than the TKE model for seabeds with submerged structures considering all the 08 cases studied. Overall performance of the TKE model is therefore found to be better than that of the Kennedy model for wave evolution over seabeds with submerged structures.

Free Water Surface Elevation
The temporal variations of free water surface are compared here across the surf zone since the wave height or mean water level only provides a phase-averaged quantity associated with wave propagation. The cases PLS03 (with plane slope) and SSB02 (with a submerged structure) are selected for discussion as they represent two different types of bottom configuration. Figures  6 and 7 depict the time series of measured and simulated free water surface elevation at nine gauge positions across the surf zone associated with case PLS03 and SSB02, respectively. The blue lines show the profiles associated with the TKE model and the red lines show those associated with the Kennedy model. Both the sub-models predict the profiles with reasonable accuracy with very little difference between them irrespective of the configuration of the bottom. However, a slight over-smoothing can be observed in profiles associated with case

Spatial and Temporal Evolution of Eddy Viscosity
The predictive skills of the two breaking models are further investigated and compared with reference to the eddy viscosity distribution across the surf zone. The cases PLS03 and SSB02 are again selected, and the temporal and spatial variations of eddy viscosity across the surf zone for each case (Figure 8) are examined. For the two cases selected, the spatial distribution of eddy viscosity at five equi-spaced phases is plotted for explaining its temporal evolution.
Similar to Shaffer et al. [24] surface roller model, the eddy viscosity production and evolution are only concentrated at the front face of the wave crest when computations are made with Kennedy model. Hence, the wave energy dissipation due to wave breaking is only concentrated at the front face of the wave crest. This is the reason why multiple peaks are observed in the eddy viscosity distribution across the surf zone. Three prominent front faces of broken wave crests can be observed approximately 2.4 m, 1.55 m and 0.8 m from the still water shoreline at the time instance T/5 for case PLS03, whereas such broken wave profiles can be observed approximately 2.3 m, 1.6 m, 0.85 m and 0.2 m from the still water shoreline at the time instant T/5 for case SSB02. The peaks in the eddy viscosity distribution follow the development and movement of surface rollers in both cases (follow 2T/5, 3T/5, 4T/5 and T). The observed wave decay is relatively smoother and gradual in PLS03 than that in SSB02 (see Figures  2 & 3). The gradual decay of wave heights is always accompanied by a gradual rise in the mean water level. On the other hand, the eddy viscosity distribution in SSB02 indicates an intense wave breaking on the submerged structure and flattening of wave profile soon after breaking, which leads to a short recovery in wave breaking. The intense breaking leads to a sudden rise in mean water level on the structure followed by a stable mean water level between the shoreline and the structure unless a secondary or multiple breaking occurs. On the contrary, the production of turbulent kinetic energy only occurs at the wave breaking point (hence, production of eddy viscosity) in TKE model. Consequently, turbulent kinetic energy is distributed across the surf zone by diffusion and advection, in addition to the relevant dissipation in the system. This leads to a very smooth eddy viscosity distribution across the surf zone. It is noted that only a single peak is observed in eddy viscosity distribution when a single breaking event occurs, but multiple peaks can be observed following multiple breaking points. In cases examined here (i.e. PLS03 and SSB02), smooth temporal and spatial distributions of eddy viscosity are observed with the peak in eddy viscosity associated with case SSB02 being steeper than that of PLS03. As explained earlier, this is due to the intense wave breaking and dissipation over the submerged structure in case SSB02. In both the cases, due to the continuous distribution of eddy viscosity across the surf zone (not just concentrated to the front face of the wave crest), an over-smoothing of wave profiles is obviously expected, i.e. the energy dissipation occurs even at the back of the wave or wave trough due to turbulence, which is not exactly false. The only concern is the magnitude of the eddy viscosity computed using the TKE model at such locations of the waves. It is noted here that model is used to compute the turbulent kinetic energy distribution in this study and separate simultaneous computations were not conducted for dissipation of the turbulent kinetic energy (i.e. − model). Although the temporal and spatial distributions of eddy viscosity in Kennedy model and TKE model differ from each other, the order of magnitude of eddy viscosity is observed to be consistent across all cases examined in this study.

Conclusions
Two eddy viscosity type wave breaking models, namely TKE model and Kennedy model (as identified in this study), investigated for their predictive skills and performances are compared by coupling them into two weakly nonlinear Boussinesq-type wave models. No effort is made to calibrate the model parameters. Therefore, Nwogu and Demirbilek [28] and Choi et al. [30] recommendations are incorporated into the models to perform all the simulations.
This comparative study demonstrates very similar results indicating relatively high accuracy in predictive skills of the two models. Kennedy model is found to be slightly superior to the TKE model when the evolution of waves over plane seabed slopes is considered. However, TKE model is found to perform slightly better than Kennedy model for wave evolution over seabed slopes with submerged structures.
Although the wave energy dissipation primarily is concentrated on the front face of the wave crests, the Kennedy model needs to incorporate possible energy dissipation over the wave profile at other locations due to the existence of turbulence. On the other hand, a mechanism must be incorporated into the turbulent kinetic energy equation to demonstrate a more realistic distribution of eddy viscosity across the surf zone. i.e. by extending the model to a superior closure model like − model.  [33], **Siddique et al. [21], ***Nwogu [25], + Hansen and Svendsen [34], ++ Ranasinghe et al. [35], +++ Ranasinghe et al. [36]