Simulation of breaking waves over different seabed configurations using phase-resolving wave models

The applicability of two phase-resolving wave models (Boussinesq-type) for simulating wave transformation over seabeds with plane slopes and with submerged structures is further investigated by coupling a turbulent kinetic energy (TKE) equation for wave breaking induced energy dissipation. Previous studies have shown that the models developed using Boussinesq equations with improved nonlinearity and frequency dispersion are capable of simulating wave evolution over mildly sloping beaches at shallow/intermediate water depths but with little emphasis on changes to bottom configurations. In this study, a published set of data is used to calibrate two phase-resolving wave models developed extending an improved form of the Boussinesq equations to include steep localised bottom slopes and wave breaking induced energy dissipation. The first model is limited to simulating wave transformation over impermeable beds but the second model is developed to incorporate porous beds. Porous damping is introduced in the second model by coupling the governing equations with a nonlinear Darcy-Forchheimer equation. The calibrated models with appropriate values for relevant model parameters are found to reproduce wave height, mean water level distributions across the surf zone for different wave conditions and bottom configurations with a high level of accuracy in 1-dimensional horizontal (1DH) wave propagation. The results obtained from this study are vital in validating a 2-dimensional horizontal (2DH) wave-current model.


INTRODUCTION
Over the years, researchers have continuously made contributions to the development and enhancement of mathematical theories for wave propagation (such as Boussinesq equations, energy balance equation, mild slope equation and RANS equations), and described wave evolutions due to various phenomena in the nearshore. Out of those main wave equations, Boussinesq-type equations (phase-resolving) have made a remarkable development over the last couple of decades largely due to the critical steps provided by Madsen et al. (1991) and Nwogu (1993). Equations by Madsen et al. (1991) were further extended in Madsen and Sorensen (1992) to include terms proportional to the bottom slope, which are essential for demonstrating shoaling characteristics of waves. Most of the subsequent theories, which were developed to extend the nonlinearity and frequency dispersion properties Madsen et al., 1996;Madsen & Schäffer,1998;Zou, 1999;Gobbi et al., 2000;Chen et al., 2003;Schäffer, 2004;Li, 2008;Karambas & Memos 2009;Chondros & Memos, 2014) were based on Madsen and Sorensen (1992) and Nwogu (1993) platform equations. Furthermore, the incorporation of porous flow resistance in Boussinesqtype models were performed by Cruz et al. (1997), Hsiao et al. (2002), Garcia et al. (2004), Chen (2006) and Cruz and Chen (2007) and this has enabled the application of such models for simulating wave propagation over permeable beds and submerged bars. However, these models are primarily limited to non-breaking wave dissipation over permeable beds; therefore, appropriate wave breaking sub-models are necessary to be incorporated to investigate wave breaking induced energy dissipation over such bottom configurations.
Modelling surf zone hydrodynamics is of great interest for many reasons. Although wave breaking is the most dominant driving mechanism for the current generation and subsequently the sediment transport in the coastal zone, due to highly nonlinear and turbulent nature of the flow, a full mathematical description of the complex hydrodynamics in the surf zone is difficult. Most of the wave theories, including Boussinesq-type, are obtained under the hypothesis of ideal and irrotational flow which restricts the direct incorporation of production and evolution of vorticity within the equations. In order to overcome this limitation, some closure models have been coupled to Boussinesq equations by researchers to simulate wave breaking and breaking induced energy dissipation. These include reasonably ad hoc additions of eddy viscosity formulations up to reasonably detailed calculations of the generation and transport of turbulent kinetic energy under the breaking wave crest or the surface rollers. 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. At a minimum, these terms are scaled similarly in order to reproduce the correct amount of energy dissipation. They are also be localised in the front face of the breaking wave, in order to provide the correct distribution of dissipation (Kirby, 2003). All these models differ on how they treat the onset and cessation of wave breaking and the rate of energy dissipation. It is also noted here that Boussinesq equations are depthintegrated equations, which limit the simulation of eddy viscosity distributions over the depth.
The first attempt to simulate wave breaking in a phase resolving model was proposed by Zelt (1991) by introducing a wave force term in the momentum equation. The term regulates the amount of energy dissipation produced by wave breaking and is related to the distribution of eddy viscosity. In the eddy viscosity formulation, to initiate and terminate the wave breaking process, some breaking detection criterion was used. Following Zelt (1991), Kennedy et al. (2000) proposed a new eddy viscosity formulation but with extensions to provide a more realistic description of the onset and cessation of breaking. However, in their formulation, no direct physical meaning can be attributed to the mixing length coefficient and quantities which govern the onset and cessation of breaking.
A simplified formulation with a more direct physical interpretation was developed and introduced into Boussinesq equations by Shaffer et al. (1993) using the surface roller concept of Svendsen (1984). They determined the variation of surface roller in space and time by adopting the heuristic geometrical approach of Deigaard (1989) with some modifications. A strong analogy can be found between Shaffer et al. (1993) and Kennedy et al. (2000) but, as it is the case with Kennedy et al. (2000), onset and cessation of wave breaking and dissipation in Shaffer et al. (1993) formulations depend on a number of parameters, which need calibration. Hence, the setting of appropriate values for these thresholds and determination of temporal and spatial evolution of surface rollers are still challenges to coastal engineering researchers. Nwogu (1996) extended the fully nonlinear Boussinesq equations of  to surf zone by introducing a term representing vertical gradient in shear stress due to breaking induced turbulent velocity fluctuations in the momentum equations. In this study, a semi-empirical turbulence closure model was used to relate the turbulence stresses to the wave orbital velocities. Following Boussinesq eddy viscosity concept, the turbulence shear stresses were approximated to be proportional to velocity gradients in vertical directions (proportionality coefficient being the eddy viscosity). This resulted in an energy dissipation term equivalent to the one proposed by Zelt (1991). However, in Nwogu (1996) formulations, the eddy viscosity is related to the turbulent kinetic energy and a turbulence length scale. Although the number of free parameters in the turbulent kinetic energy equation is relatively less, there is still an ambiguity in defining an appropriate value for the turbulence length scale (mixing length).
A few studies can be found in the literature on modelling of breaking waves over permeable bottoms/ structures using Boussinesq-type models. Johnson et al. (2005) used a higher-order Boussinesq-type model in 2DH while Metallinos et al. (2016) used the fully dispersive and highly nonlinear Boussinesq-type model of Chondros and Memos (2014) to simulate the wave propagation over a submerged permeable structure in 1DH. In the presence of permeable structures, both models incorporated two additional terms accounting for the interaction between the wave motion outside and the flow within the rubble mound, one in the continuity equation and one in the momentum equation in accordance to the model's extension provided by Cruz et al. (1997). Both models used lower-order porous damping terms and eddy viscosity production and wave breaking induced energy dissipation simulations were performed using formulations proposed by Kennedy et al. (2000).
Journal of the National Science Foundation of Sri Lanka 49 (1) March 2021 In the present study, the wave transformation in nearshore is modelled using a truncated form of Chen (2006) and Nwogu (1993) Boussinesq-type wave equations coupled with one-equation turbulence model (Nowgu, 1996) type) to determine the temporal and spatial evolution of the turbulent kinetic energy produced by wave breaking. Nwogu (1996), Tajima et al. (2007), Ranasinghe et al. (2009aRanasinghe et al. ( , 2011, Siddique et al. (2017) and Ranasinghe (2018) followed a similar approach for wave breaking induced energy dissipation over plane seabed slopes and seabeds with submerged structures; however, all used values for mixing length which are specific to the cases they studied. BOUSS-2D, a Boussinesq wave model developed to simulate waves and currents in nearshore regions and harbours by US Army Corps of Engineers (Nwogu & Demirbilek, 2001) too employs a TKE equation for wave breaking induced energy dissipation and recommends a value equal to water wave height for mixing length but with no emphasis on any variation due to bottom configuration or wave period. Hence, considering the ambiguities in setting appropriate values for wave breaking parameters for this type of numerical model, an attempt is made in this study to revisit the incorporation of TKE equation for eddy viscosity formulation, set a simple wave breaking detection criterion for two weakly nonlinear Boussinesqtype models, and propose an appropriate representation for mixing length using a published set of data covering plane seabed slopes and seabeds with submerged structures for 1-dimensional horizontal (1DH) wave propagation. These results are expected to be used to validate a 2-dimensional horizontal (2DH) Boussinesqtype wave-current model.

The governing equations
Although a number of upgrades have been developed over the past two decades, still Nwogu (1993) and Madsen and Sorensen (1992) Boussineq-type equations are preferred and used for commercial purposes in nearshore wave current modelling due to simplicity and stability of these models when used to simulate waves and currents over mildly sloping impervious beaches. The models developed with fully nonlinear dispersive Boussinesq-type equations or Navier-Stokes equations are difficult to be deployed in practical applications due to stability issues and very high computational cost.
As highlighted earlier, a truncated form of Chen (2006) and Nwogu (1993) Boussinesq-type equations are used for the investigations in the present study. The truncated Chen (2006) equations are analogous to Nwogu (1993) equations when porous damping terms are switched off. A short description of the governing equations and different sub-models used in the two models is presented below.
In one horizontal dimension (1DH), the truncated form of Chen (2006) equations (i.e. by dropping higherorder nonlinear and porous damping terms) read as:

Continuity equation:
where, = ℎ + + 1 ℎ 2 + 2 ℎ ℎ + ℎ = ℎ + ℎ 2 3 ℎ − 4 ℎ + 4 ℎ 2 ℎ 1 = where is the porosity of the porous material, is the free water s water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total particle velocity at reference water depth , and is the water part depth inside porous layer along the wave propagation direction (re and are determined by matching dispersion characteristics of wave porous layer (not discussed in this paper).
is where is the porosity of the porous material, is the free water surface e water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total water particle velocity at reference water depth , and is the water particle velo depth inside porous layer along the wave propagation direction (refer to Fi and are determined by matching dispersion characteristics of wave motion i porous layer (not discussed in this paper).
is the porous resistance force follows; = + + where and are respectively, the linear and nonlinear porous resistance the added mass coefficient. The porous resistance coefficients are estimate relationships following Sollitt & Cross (1972): ... where is the porosity of the porous material, is the free water surface e water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total water particle velocity at reference water depth , and is the water particle velo depth inside porous layer along the wave propagation direction (refer to Fi and are determined by matching dispersion characteristics of wave motion i porous layer (not discussed in this paper).
is where is the porosity of the porous material, is the free water surface water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total wate particle velocity at reference water depth , and is the water particle ve depth inside porous layer along the wave propagation direction (refer to and are determined by matching dispersion characteristics of wave motion porous layer (not discussed in this paper).
is the porous resistance forc follows; where, = ℎ + + 1 ℎ 2 + 2 ℎ ℎ + ℎ = ℎ + ℎ 2 3 ℎ − 4 ℎ + 4 ℎ 2 ℎ 1 = where is the porosity of the porous material, is the free water surface elevat water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total water depth particle velocity at reference water depth , and is the water particle velocity a depth inside porous layer along the wave propagation direction (refer to Figure   and are determined by matching dispersion characteristics of wave motion in bot porous layer (not discussed in this paper).
is where is the porosity of the porous material, is the free water surface elev water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total water dep particle velocity at reference water depth , and is the water particle velocit depth inside porous layer along the wave propagation direction (refer to Figur and are determined by matching dispersion characteristics of wave motion in b porous layer (not discussed in this paper).
is the porous resistance force eva follows; where, = ℎ + + 1 ℎ 2 + 2 ℎ ℎ + ℎ = ℎ + ℎ 2 3 ℎ − 4 ℎ + 4 ℎ 2 ℎ 1 = where is the porosity of the porous material, is the free water surface elevation, water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total water depth, particle velocity at reference water depth , and is the water particle velocity at re depth inside porous layer along the wave propagation direction (refer to Figure 1). and are determined by matching dispersion characteristics of wave motion in both c porous layer (not discussed in this paper).
is where is the porosity of the porous material, is the free water s water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total particle velocity at reference water depth , and is the water part depth inside porous layer along the wave propagation direction (re and are determined by matching dispersion characteristics of wave porous layer (not discussed in this paper).
is the porous resistanc follows; where, = ℎ + + 1 ℎ 2 + 2 ℎ ℎ + ℎ = ℎ + ℎ 2 3 ℎ − 4 ℎ + 4 ℎ 2 ℎ 1 = where is the porosity of the porous material, is the free water surface ele water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total water d particle velocity at reference water depth , and is the water particle veloc depth inside porous layer along the wave propagation direction (refer to is the porous resistance force e follows; ... (5) Momentum equation for clear water: where, = ℎ + + 1 ℎ 2 + 2 ℎ ℎ + ℎ = ℎ + ℎ 2 3 ℎ − 4 ℎ + 4 ℎ 2 ℎ 1 = where is the porosity of the porous material, is the free water s water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total particle velocity at reference water depth , and is the water part depth inside porous layer along the wave propagation direction (re and are determined by matching dispersion characteristics of wave porous layer (not discussed in this paper).
is the porous resistanc follows; where, Momentum equation for clear water: Momentum equation for porous layer: (8) 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 (refer to Figure 1). where is the porosity of the porous material, is the free water s water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total particle velocity at reference water depth , and is the water part depth inside porous layer along the wave propagation direction (re and are determined by matching dispersion characteristics of wave porous layer (not discussed in this paper).
is the porous resistanc ... (7) Momentum equation for porous layer: where, = ℎ + + 1 ℎ 2 + 2 ℎ ℎ + ℎ = ℎ + ℎ 2 3 ℎ − 4 ℎ + 4 ℎ 2 ℎ 1 = where is the porosity of the porous material, is the free water surface water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total water particle velocity at reference water depth , and is the water particle vel depth inside porous layer along the wave propagation direction (refer to F and are determined by matching dispersion characteristics of wave motion porous layer (not discussed in this paper).
is the porous resistance force follows; where, = ℎ + + 1 ℎ 2 + 2 ℎ ℎ + ℎ = ℎ + ℎ 2 3 ℎ − 4 ℎ + 4 ℎ 2 ℎ 1 = where is the porosity of the porous material, is the free water surface water depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total water particle velocity at reference water depth , and is the water particle vel depth inside porous layer along the wave propagation direction (refer to F and are determined by matching dispersion characteristics of wave motion porous layer (not discussed in this paper).
is the porous resistance force follows; ... (8)  (8) 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 (refer to Figure 1). 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 (1972): is the porosity of the porous material, 5 + + + 2 2 + ℎ + 1 2 ℎ 2 − ℎ ℎ − ℎ ℎ − 1 2 ℎ 2 + ℎ ℎ + + 2 2 + ℎ + ℎ 2 2 + ℎ ℎ = 0 (8) 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 (refer to Figure 1). 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 (1972): is the free water surface elevation, is the free water surface elevation, ℎ is the clear ℎ ( ℎ + ℎ ) is the total water depth, is the water , and is the water particle velocity at reference water propagation direction (refer to Figure 1). The factors, n characteristics of wave motion in both clear water and is the porous resistance force evaluated at = as r and nonlinear porous resistance coefficients, and is sistance coefficients are estimated from the following 2): is the clear water depth, 5 + + + 2 2 + ℎ + 1 2 ℎ 2 − ℎ ℎ − ℎ ℎ − 1 2 ℎ 2 + ℎ ℎ + + 2 2 + ℎ + ℎ 2 2 + ℎ ℎ = 0 (8) 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 (refer to Figure 1). 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 (1972): is the porous layer depth, 5 + + + 2 2 + ℎ + 1 2 ℎ 2 − ℎ ℎ − ℎ ℎ − 1 2 ℎ 2 + ℎ ℎ + + 2 2 + ℎ + ℎ 2 2 + ℎ ℎ = 0 (8) 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 (refer to Figure 1). 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 (1972): is the total water depth − ℎ ℎ − ℎ ℎ − + ℎ 2 2 + ℎ ℎ = 0 (8) free water surface elevation, ℎ is the clear is the total water depth, is the water e water particle velocity at reference water direction (refer to Figure 1). The factors, ics of wave motion in both clear water and us resistance force evaluated at = as (9) r porous resistance coefficients, and is ficients are estimated from the following is the water particle velocity at reference water depth 5 uation for porous layer: (8) porosity of the porous material, is the free water surface elevation, ℎ is the clear is the porous layer depth, ℎ ( ℎ + ℎ ) is the total water depth, is the water ty at reference water depth , and is the water particle velocity at reference water e porous layer along the wave propagation direction (refer to Figure 1). The factors, rmined by matching dispersion characteristics of wave motion in both clear water and not discussed in this paper). is the porous resistance force evaluated at = as are respectively, the linear and nonlinear porous resistance coefficients, and is ss coefficient. The porous resistance coefficients are estimated from the following ollowing Sollitt & Cross (1972): , and 5 uation for porous layer: (8) porosity of the porous material, is the free water surface elevation, ℎ is the clear is the porous layer depth, ℎ ( ℎ + ℎ ) is the total water depth, is the water ty at reference water depth , and is the water particle velocity at reference water porous layer along the wave propagation direction (refer to Figure 1). The factors, rmined by matching dispersion characteristics of wave motion in both clear water and not discussed in this paper). is the porous resistance force evaluated at = as are respectively, the linear and nonlinear porous resistance coefficients, and is ss coefficient. The porous resistance coefficients are estimated from the following ollowing Sollitt & Cross (1972): is the water particle velocity at reference water depth 5 + + + 2 2 + ℎ + 1 2 ℎ 2 − ℎ ℎ − ℎ ℎ − 1 2 ℎ 2 + ℎ ℎ + + 2 2 + ℎ + ℎ 2 2 + ℎ ℎ = 0 (8) 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 (refer to Figure 1). 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 (1972): inside porous layer along the wave propagation direction (Figure 1). The factors, 5 er: is the free water surface elevation, ℎ is the clear er depth, ℎ ( ℎ + ℎ ) is the total water depth, is the water depth , and is the water particle velocity at reference water the wave propagation direction (refer to Figure 1). The factors, dispersion characteristics of wave motion in both clear water and paper). is the porous resistance force evaluated at = as the linear and nonlinear porous resistance coefficients, and is orous resistance coefficients are estimated from the following ross (1972): (8) 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 (refer to Figure 1). 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 (1972): are determined by matching dispersion characteristics of wave motion in both clear water and porous layer (not discussed in this paper). (8) e is the porosity of the porous material, is the free water surface elevation, ℎ is the clear r depth, ℎ is the porous layer depth, ℎ ( ℎ + ℎ ) is the total water depth, is the water le velocity at reference water depth , and is the water particle velocity at reference water inside porous layer along the wave propagation direction (refer to Figure 1). The factors, are determined by matching dispersion characteristics of wave motion in both clear water and s layer (not discussed in this paper). is the porous resistance force evaluated at = as ws; e and are respectively, the linear and nonlinear porous resistance coefficients, and is dded mass coefficient. The porous resistance coefficients are estimated from the following onships following Sollitt & Cross (1972): is the porous resistance force evaluated at ℎ − ℎ ℎ − ℎ 2 2 + ℎ ℎ = 0 (8) ater surface elevation, ℎ is the clear e total water depth, is the water er particle velocity at reference water ion (refer to Figure 1). The factors, wave motion in both clear water and sistance force evaluated at = as (9) ous resistance coefficients, and is ts are estimated from the following as follows; (8) 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 (refer to Figure 1). 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 (1972): ... (9) where 5 + + + 2 2 + ℎ + 1 2 ℎ 2 − ℎ ℎ − ℎ ℎ − 1 2 ℎ 2 + ℎ ℎ + + 2 2 + ℎ + ℎ 2 2 + ℎ ℎ = 0 (8) 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 (refer to Figure 1). 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 (1972): and 5 + + + 2 2 + ℎ + 1 2 ℎ 2 − ℎ ℎ − ℎ ℎ − 1 2 ℎ 2 + ℎ ℎ + + 2 2 + ℎ + ℎ 2 2 + ℎ ℎ = 0 (8) 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 (refer to Figure 1). 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 (1972): are respectively, the linear and nonlinear porous resistance coefficients, and 5 er: is the free water surface elevation, ℎ is the clear er depth, ℎ ( ℎ + ℎ ) is the total water depth, is the water depth , and is the water particle velocity at reference water the wave propagation direction (refer to Figure 1). The factors, dispersion characteristics of wave motion in both clear water and paper). is the porous resistance force evaluated at = as (9) the linear and nonlinear porous resistance coefficients, and is orous resistance coefficients are estimated from the following ross (1972): is the added mass coefficient. The porous resistance coefficients are estimated from the following relationships following Sollitt & Cross (1972) 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 (1995).
where 50 is the median diameter of the porous material, , and are empirical coefficients that, in principle, are determined by experiments. However, van Gent (1995) 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 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. where is the kinematic viscosity of water (∽ 10 −6 m 2 /s). The intri mass coefficient, , and nonlinear drag coefficient, , are determ proposed by van Gent (1995).
where 50 is the median diameter of the porous material, , and that, in principle, are determined by experiments. However, van G although the coefficients and depend on the parameters lik orientation of the porous material, the values 1000 and 1.1 can be characteristic length scale 50 is used to define the size of poro suggested using a constant value for empirical coefficient, of 0.34 complex dependency of coefficient, on the flow field. The poros used in the present study was found to be 0.44. The above equations yield to Nwogu (1993) Boussinesq-type equatio terms are ignored. i.e., Continuity equation: Momentum equation for clear water: where is the kinematic viscosity of water (∽ 10 −6 m 2 /s). The intrinsic mass coefficient, , and nonlinear drag coefficient, , are determine proposed by van Gent (1995).
where 50 is the median diameter of the porous material, , and that, in principle, are determined by experiments. However, van Gen although the coefficients and depend on the parameters like s orientation of the porous material, the values 1000 and 1.1 can be use characteristic length scale 50 is used to define the size of porous suggested using a constant value for empirical coefficient, of 0.34 for complex dependency of coefficient, on the flow field. The porosity, used in the present study was found to be 0.44. The above equations yield to Nwogu (1993) Boussinesq-type equations w terms are ignored. i.e., Continuity equation: Momentum equation for clear water: where is the kinematic viscosity of water (∽ 10 −6 m 2 /s). The intrinsic pe mass coefficient, , and nonlinear drag coefficient, , are determined b proposed by van Gent (1995).
where 50 is the median diameter of the porous material, , and are that, in principle, are determined by experiments. However, van Gent (1 although the coefficients and depend on the parameters like shap orientation of the porous material, the values 1000 and 1.1 can be used f characteristic length scale 50 is used to define the size of porous ma suggested using a constant value for empirical coefficient, of 0.34 for sim complex dependency of coefficient, on the flow field. The porosity, n o used in the present study was found to be 0.44. The above equations yield to Nwogu (1993) Boussinesq-type equations whe terms are ignored. i.e., Continuity equation: Momentum equation for clear water: ... (11) where 6 where is the kinematic viscosity of water (∽ 10 −6 m 2 /s). The mass coefficient, , and nonlinear drag coefficient, , are de proposed by van Gent (1995).
where 50 is the median diameter of the porous material, , a that, in principle, are determined by experiments. However, va although the coefficients and depend on the parameters orientation of the porous material, the values 1000 and 1.1 can characteristic length scale 50 is used to define the size of p suggested using a constant value for empirical coefficient, of 0 complex dependency of coefficient, on the flow field. The p used in the present study was found to be 0.44. The above equations yield to Nwogu (1993) Boussinesq-type equ terms are ignored. i.e., Continuity equation: Momentum equation for clear water: is the median diameter of the porous material, 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 (1995).
where 50 is the median diameter of the porous material, , and are empirical coefficients that, in principle, are determined by experiments. However, van Gent (1995) 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 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 (1993) Boussinesq-type equations when the porous damping terms are ignored. i.e., Continuity equation: Momentum equation for clear water: 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 (1995).
where 50 is the median diameter of the porous material, , and are empirical coefficients that, in principle, are determined by experiments. However, van Gent (1995) 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 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 (1993) Boussinesq-type equations when the porous damping terms are ignored. i.e., Continuity equation: Momentum equation for clear water: are empirical coefficients that, in principle, are determined by experiments. However, van Gent (1995) concluded that, although the coefficients where is the kinematic viscosity of water (∽ 10 −6 m 2 /s). The intrinsic permeability, , add mass coefficient, , and nonlinear drag coefficient, , are determined by empirical formu proposed by van Gent (1995).
where 50 is the median diameter of the porous material, , and are empirical coefficie that, in principle, are determined by experiments. However, van Gent (1995) concluded th although the coefficients and depend on the parameters like shape, aspect ratio or orientation of the porous material, the values 1000 and 1.1 can be used for respectively, if characteristic length scale 50 is used to define the size of porous material. In addition, suggested using a constant value for empirical coefficient, of 0.34 for simplicity neglecting complex dependency of coefficient, on the flow field. The porosity, n of the porous mater used in the present study was found to be 0.44. The above equations yield to Nwogu (1993) Boussinesq-type equations when the porous dampi terms are ignored. i.e., Continuity equation: Momentum equation for clear water: where is the kinematic viscosity of water (∽ 10 −6 m 2 /s). The intrinsic permeability, , add mass coefficient, , and nonlinear drag coefficient, , are determined by empirical formul proposed by van Gent (1995).
where 50 is the median diameter of the porous material, , and are empirical coefficien that, in principle, are determined by experiments. However, van Gent (1995) concluded th although the coefficients and depend on the parameters like shape, aspect ratio or t orientation of the porous material, the values 1000 and 1.1 can be used for respectively, if t characteristic length scale 50 is used to define the size of porous material. In addition, suggested using a constant value for empirical coefficient, of 0.34 for simplicity neglecting t complex dependency of coefficient, on the flow field. The porosity, n of the porous materi used in the present study was found to be 0.44. The above equations yield to Nwogu (1993) Boussinesq-type equations when the porous dampin terms are ignored. i.e., Continuity equation: Momentum equation for clear water: depend on the parameters like shape, aspect ratio or the orientation of the porous material, the values 1000 and 1.1 can be used respectively, if the characteristic length scale where is the kinematic viscosity of water (∽ 10 −6 m 2 /s). The intrinsic permeabil mass coefficient, , and nonlinear drag coefficient, , are determined by empi proposed by van Gent (1995).
where 50 is the median diameter of the porous material, , and are empiric that, in principle, are determined by experiments. However, van Gent (1995) c although the coefficients and depend on the parameters like shape, aspec orientation of the porous material, the values 1000 and 1.1 can be used for respe characteristic length scale 50 is used to define the size of porous material. In suggested using a constant value for empirical coefficient, of 0.34 for simplicity complex dependency of coefficient, on the flow field. The porosity, n of the p used in the present study was found to be 0.44. The above equations yield to Nwogu (1993) Boussinesq-type equations when the po terms are ignored. i.e., Continuity equation: Momentum equation for clear water: is used to define the size of porous material. In addition, he suggested using a constant value (0.34) for empirical coefficient, where is the kinematic viscosity of water (∽ 10 −6 m 2 /s). The intrinsic permeability, mass coefficient, , and nonlinear drag coefficient, , are determined by empirical fo proposed by van Gent (1995).
where 50 is the median diameter of the porous material, , and are empirical coef that, in principle, are determined by experiments. However, van Gent (1995) conclude although the coefficients and depend on the parameters like shape, aspect ratio orientation of the porous material, the values 1000 and 1.1 can be used for respectively characteristic length scale 50 is used to define the size of porous material. In addit suggested using a constant value for empirical coefficient, of 0.34 for simplicity neglec complex dependency of coefficient, on the flow field. The porosity, n of the porous m used in the present study was found to be 0.44.  (Ranasinghe et al., 2009b) The above equations yield to Nwogu (1993) Boussinesq-type equations when the porous d terms are ignored. i.e., Continuity equation: Momentum equation for clear water: for simplicity neglecting the complex dependency of coefficient, where is the kinematic viscosity of water (∽ 10 −6 m 2 /s mass coefficient, , and nonlinear drag coefficient, proposed by van Gent (1995).
where 50 is the median diameter of the porous material, that, in principle, are determined by experiments. How although the coefficients and depend on the par orientation of the porous material, the values 1000 and characteristic length scale 50 is used to define the s suggested using a constant value for empirical coefficient, complex dependency of coefficient, on the flow field used in the present study was found to be 0.44. The above equations yield to Nwogu (1993) Boussinesq-t terms are ignored. i.e., Continuity equation: Momentum equation for clear water: 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 (1993) Boussinesq-type equations when the porous damping terms are ignored. i.e.,

Continuity equation:
6 that, in principle, are determined by experiments. However, van although the coefficients and depend on the parameters l orientation of the porous material, the values 1000 and 1.1 can be characteristic length scale 50 is used to define the size of por suggested using a constant value for empirical coefficient, of 0.3 complex dependency of coefficient, on the flow field. The poro used in the present study was found to be 0.44. The above equations yield to Nwogu (1993) Boussinesq-type equati terms are ignored. i.e., Continuity equation: Momentum equation for clear water: ... (12) Momentum equation for clear water: Energy Dissipation due to Friction and Wave Breaking Two additional terms are introduced into momentum conservatio (2006) type is used and Equation 13 when Nwogu (1993) type of wave energy due to bottom friction and wave breaking, which = 1 ℎ+ where is the wave friction factor, and is the near-bottom axis. Following Kennedy et al. (2000), a simple eddy viscosity simulate energy dissipation due to wave breaking induced ene momentum-mixing term into momentum conservation equation The rate of wave energy dissipation is expected to be govern viscosity, , which is related to the turbulent kinetic energy, The turbulent kinetic energy is determined from a semi-empiric term for turbulent kinetic energy production by wave breaking, i model (Nwogu (1996)): ...(13)

Energy dissipation due to friction and wave breaking
Two additional terms are introduced into momentum conservation equation [Equation 6 when Chen (2006) type is used and Equation 13 when Nwogu (1993) type is used] to simulate the dissipation of wave energy due to bottom friction and wave breaking, which are represented by is the wave friction factor, and is the near-bottom wave orbital velocity along the xollowing Kennedy et al. (2000), a simple eddy viscosity type of formulation is adopted to te energy dissipation due to wave breaking induced energy dissipation (by introducing the ntum-mixing term into momentum conservation equation in the x-direction).
te of wave energy dissipation is expected to be governed by the magnitude of the eddy ity, , which is related to the turbulent kinetic energy, , and a turbulent length scale, . rbulent kinetic energy is determined from a semi-empirical transport equation with a source , respectively.
Energy Dissipation due to Friction and Wave Breaking Two additional terms are introduced into momentum conservation equ (2006) type is used and Equation 13 when Nwogu (1993) type is used of wave energy due to bottom friction and wave breaking, which are re = 1 ℎ+ where is the wave friction factor, and is the near-bottom wave axis. Following Kennedy et al. (2000), a simple eddy viscosity type simulate energy dissipation due to wave breaking induced energy dis momentum-mixing term into momentum conservation equation in the x = 1 ℎ+ ℎ + The rate of wave energy dissipation is expected to be governed by viscosity, , which is related to the turbulent kinetic energy, , and The turbulent kinetic energy is determined from a semi-empirical tran term for turbulent kinetic energy production by wave breaking, i.e., on ... (14) where + + + 1 ℎ 2 + 2 ℎ ℎ = 0 Energy Dissipation due to Friction and Wave Breaking Two additional terms are introduced into momentum conservation e (2006) type is used and Equation 13 when Nwogu (1993) type is u of wave energy due to bottom friction and wave breaking, which are = 1 ℎ+ where is the wave friction factor, and is the near-bottom wa axis. Following Kennedy et al. (2000), a simple eddy viscosity ty simulate energy dissipation due to wave breaking induced energy momentum-mixing term into momentum conservation equation in t = 1 ℎ+ ℎ + The rate of wave energy dissipation is expected to be governed viscosity, , which is related to the turbulent kinetic energy, , a The turbulent kinetic energy is determined from a semi-empirical is the wave friction factor, and Energy Dissipation due to Friction and Wave Breaking Two additional terms are introduced into momentum conservation e (2006) type is used and Equation 13 when Nwogu (1993) type is of wave energy due to bottom friction and wave breaking, which are = 1 ℎ+ where is the wave friction factor, and is the near-bottom w axis. Following Kennedy et al. (2000), a simple eddy viscosity ty simulate energy dissipation due to wave breaking induced energy momentum-mixing term into momentum conservation equation in t = 1 ℎ+ ℎ + The rate of wave energy dissipation is expected to be governed viscosity, , which is related to the turbulent kinetic energy, , a The turbulent kinetic energy is determined from a semi-empirical is the near-bottom wave orbital velocity along the x-axis. Following Kennedy et al. (2000), a simple eddy 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 (1995).
where 50 is the median diameter of the porous material, , and are empirical coefficients that, in principle, are determined by experiments. However, van Gent (1995) 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 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.  (Ranasinghe et al., 2009b) The above equations yield to Nwogu (1993) Boussinesq-type equations when the porous damping terms are ignored. i.e., Continuity equation: Momentum equation for clear water: 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 (1995).
where 50 is the median diameter of the porous material, , and are empirical coefficients that, in principle, are determined by experiments. However, van Gent (1995) 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 respectively, if the characteristic length scale 50 is used to define the size of porous material. In addition, he is the kinematic viscosity of water (10) osity of water (∽ 10 −6 m 2 /s). The intrinsic permeability, , added nlinear drag coefficient, , are determined by empirical formulae eter of the porous material, , and are empirical coefficients ined by experiments. However, van Gent (1995) concluded that, and depend on the parameters like shape, aspect ratio or the terial, the values 1000 and 1.1 can be used for respectively, if the is used to define the size of porous material. In addition, he . The intrinsic permeability, osity of water (∽ 10 −6 m 2 /s). The intrinsic permeability, , added onlinear drag coefficient, , are determined by empirical formulae eter of the porous material, , and are empirical coefficients ined by experiments. However, van Gent (1995) concluded that, and depend on the parameters like shape, aspect ratio or the terial, the values 1000 and 1.1 can be used for respectively, if the 50 is used to define the size of porous material. In addition, he added mass coefficient, 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 (1995).
where 50 is the median diameter of the porous material, , and are empirical coefficients that, in principle, are determined by experiments. However, van Gent (1995) 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 respectively, if the characteristic length scale 50 is used to define the size of porous material. In addition, he and nonlinear drag coefficient, 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 (1995).
where 50 is the median diameter of the porous material, , and are empirical coefficients that, in principle, are determined by experiments. However, van Gent (1995) 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 respectively, if the characteristic length scale 50 is used to define the size of porous material. In addition, he , are determined by empirical formulae proposed by van Gent (1995).
Journal of the National Science Foundation of Sri Lanka 49 (1) March 2021 viscosity type of formulation is adopted to simulate energy dissipation due to wave breaking induced energy dissipation (by introducing the momentum-mixing term into momentum conservation equation in the x-direction).
where is the wave friction factor, and is the near-bottom wave orbital velocity along the xaxis. Following Kennedy et al. (2000), a simple eddy viscosity type of formulation is adopted to simulate energy dissipation due to wave breaking induced energy dissipation (by introducing the momentum-mixing term into momentum conservation equation in the x-direction).
The rate of wave energy dissipation is expected to be governed by the magnitude of the eddy viscosity, , which is related to the turbulent kinetic energy, , 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 (Nwogu (1996)): 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 horizontal velocity at a reference water depth, , exceeds a threshold velocity (a fraction of the wave celerity, C). The turbulent length scale, is to be quantified appropriately as part of the study, while 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, and that is obtained by assuming a local balance between production and dissipation of turbulent kinetic energy as follows: ... (15) The rate of wave energy dissipation is expected to be governed by the magnitude of the eddy viscosity, v e , which is related to the turbulent kinetic energy, k, and a turbulent length scale, l m . 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 (Nwogu (1996): 7 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. (2000), a simple eddy viscosity type of formulation is adopted to simulate energy dissipation due to wave breaking induced energy dissipation (by introducing the momentum-mixing term into momentum conservation equation in the x-direction).
The rate of wave energy dissipation is expected to be governed by the magnitude of the eddy viscosity, , which is related to the turbulent kinetic energy, , 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 (Nwogu (1996)): 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 horizontal velocity at a reference water depth, , exceeds a threshold velocity (a fraction of the wave celerity, C). The turbulent length scale, is to be quantified appropriately as part of the study, while 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, and that is obtained by assuming a local balance between production and dissipation of turbulent kinetic energy as follows:  Nwogu (1993) 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. (2000), a simple eddy viscosity type of formulation is adopted to simulate energy dissipation due to wave breaking induced energy dissipation (by introducing the momentum-mixing term into momentum conservation equation in the x-direction).
The rate of wave energy dissipation is expected to be governed by the magnitude of the eddy viscosity, , which is related to the turbulent kinetic energy, , 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 (Nwogu (1996)): 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 horizontal velocity at a reference water depth, , exceeds a threshold velocity (a fraction of the wave celerity, C). The turbulent length scale, is to be quantified appropriately as part of the study, while 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, and that is obtained by assuming a local balance between production and dissipation of turbulent kinetic energy as follows: ... (17) where 7 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. (2000), a simple eddy viscosity type of formulation is adopted to simulate energy dissipation due to wave breaking induced energy dissipation (by introducing the momentum-mixing term into momentum conservation equation in the x-direction).
The rate of wave energy dissipation is expected to be governed by the magnitude of the eddy viscosity, , which is related to the turbulent kinetic energy, , 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 (Nwogu (1996)): 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 horizontal velocity at a reference water depth, , exceeds a threshold velocity (a fraction of the wave celerity, C). The turbulent length scale, is to be quantified appropriately as part of the study, while 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, and that is obtained by assuming a local balance between production and dissipation of turbulent kinetic energy as follows: is the free surface wave orbital velocity along the x-axis. The parameter is introduced to ensure that the turbulence is produced only when horizontal velocity at a reference water depth, 7 (2006) type is used and Equation 13 when Nwogu (1993) type is used] to simulate the d 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 alo axis. Following Kennedy et al. (2000), a simple eddy viscosity type of formulation is ad simulate energy dissipation due to wave breaking induced energy dissipation (by introd momentum-mixing term into momentum conservation equation in the x-direction).
The rate of wave energy dissipation is expected to be governed by the magnitude of viscosity, , which is related to the turbulent kinetic energy, , and a turbulent length s The turbulent kinetic energy is determined from a semi-empirical transport equation with term for turbulent kinetic energy production by wave breaking, i.e., one equation turbulenc model (Nwogu (1996)): where is the free surface wave orbital velocity along the x-axis. The parameter is in to ensure that the turbulence is produced only when horizontal velocity at a reference wa , exceeds a threshold velocity (a fraction of the wave celerity, C). The turbulent leng is to be quantified appropriately as part of the study, while the empirical coefficients, for this sub-model are set to be 0.08 and 1 respectively. To trigger the turbulent kinet equation, an initial estimate of eddy viscosity, is necessary, and that is obtained by as local balance between production and dissipation of turbulent kinetic energy as follows:  Nwogu (1993) type is used] to simulate the dissipatio 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 axis. Following Kennedy et al. (2000), a simple eddy viscosity type of formulation is adopted simulate energy dissipation due to wave breaking induced energy dissipation (by introducing th momentum-mixing term into momentum conservation equation in the x-direction).
The rate of wave energy dissipation is expected to be governed by the magnitude of the edd viscosity, , which is related to the turbulent kinetic energy, , and a turbulent length scale, The turbulent kinetic energy is determined from a semi-empirical transport equation with a sour term for turbulent kinetic energy production by wave breaking, i.e., one equation turbulence closu model (Nwogu (1996)): where is the free surface wave orbital velocity along the x-axis. The parameter is introduce to ensure that the turbulence is produced only when horizontal velocity at a reference water dept , exceeds a threshold velocity (a fraction of the wave celerity, C). The turbulent length scal is to be quantified appropriately as part of the study, while the empirical coefficients, , and for this sub-model are set to be 0.08 and 1 respectively. To trigger the turbulent kinetic energ equation, an initial estimate of eddy viscosity, is necessary, and that is obtained by assuming local balance between production and dissipation of turbulent kinetic energy as follows: is to be quantified appropriately as part of the study, while the empirical coefficients, tion equation [Equation 6 when Chen e is used] to simulate the dissipation ch are represented by and .
m wave orbital velocity along the xity type of formulation is adopted to nergy dissipation (by introducing the n in the x-direction).
-axis. The parameter is introduced al velocity at a reference water depth, erity, C). The turbulent length scale, the empirical coefficients, , and trigger the turbulent kinetic energy y, and that is obtained by assuming a kinetic energy as follows: (18) 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, 7 Energy Dissipation due to Friction and Wave Breaking Two additional terms are introduced into momentum conservation equation [Equation 6 when Chen (2006) type is used and Equation 13 when Nwogu (1993) 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. (2000), a simple eddy viscosity type of formulation is adopted to simulate energy dissipation due to wave breaking induced energy dissipation (by introducing the momentum-mixing term into momentum conservation equation in the x-direction).
The rate of wave energy dissipation is expected to be governed by the magnitude of the eddy viscosity, , which is related to the turbulent kinetic energy, , 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 (Nwogu (1996)): 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 horizontal velocity at a reference water depth, , exceeds a threshold velocity (a fraction of the wave celerity, C). The turbulent length scale, is to be quantified appropriately as part of the study, while 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, and that is obtained by assuming a local balance between production and dissipation of turbulent kinetic energy as follows: is necessary, and that is obtained by assuming a local balance between production and dissipation of turbulent kinetic energy as follows: 7 Energy Dissipation due to Friction and Wave Breaking Two additional terms are introduced into momentum conservation equation [Equation 6 when Chen (2006) type is used and Equation 13 when Nwogu (1993) type is used] to simulate the dissipation of wave energy due to bottom friction and wave breaking, which are represented by and .
= 1 ℎ+ (14) where is the wave friction factor, and is the near-bottom wave orbital velocity along the xaxis. Following Kennedy et al. (2000), a simple eddy viscosity type of formulation is adopted to simulate energy dissipation due to wave breaking induced energy dissipation (by introducing the momentum-mixing term into momentum conservation equation in the x-direction).
The rate of wave energy dissipation is expected to be governed by the magnitude of the eddy viscosity, , which is related to the turbulent kinetic energy, , 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 (Nwogu (1996)): 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 horizontal velocity at a reference water depth, , exceeds a threshold velocity (a fraction of the wave celerity, C). The turbulent length scale, is to be quantified appropriately as part of the study, while 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, and that is obtained by assuming a local balance between production and dissipation of turbulent kinetic energy as follows:

Offshore open boundary
As there is significant wave reflection from submerged breakwaters, it is essential to absorb the waves, which propagate back to wave incident boundary to prevent a build-up of wave energy inside the computational domain. This is achieved with the line boundary method proposed by Ishii et al. (1994). In this method, the neighbouring grid points are divided into two types according to their position with respect to line boundary. Both incident and outgoing waves exist inside line boundary (region 1), whereas only outgoing waves exist outside the line boundary (region 2) as shown in Figure 2. To completely dissipate the energy of outgoing porosity, possible wetting and drying coexisting field (i.e. the wave trough level was allowed to touch the crown of the structure) was allowed in simulations associated with solid submerged structures following Kennedy et al. (2000). However, this was not permitted in simulations associated with porous submerged structures due to the clear existence of water and porous layers in the model.

Numerical scheme
The truncated Chen (2006) and Nwogu (1993) 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 discretised with grid size, in the -direction and the equation variables and are defined at grid points in a staggered scheme as shown in Figure 2. The water depth and surface elevation are defined at grid points , while velocities are defined half a grid point on either side of elevation grid points. The external boundaries of the computational domain correspond to velocity grid points. The first-order spatial derivatives are differenced to , which automatically eliminates error terms that would be of the same form as dispersive terms ]. The secondorder spatial derivatives are discretized to discretized to Δ 4 , while the advectio difference scheme. The computations are carr for clear water, , reference velocity for simultaneously.

Experimental Set-up
The numerical models were first validated an data associated with 15 cases. They inclu submerged structures and porous submerged models were constructed with concrete and , while the advection terms are differenced with second-order upwind difference scheme. The computations are carried out for free surface elevation, 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, , reference velocity for porous medium, and turbulent kinetic energy, k simultaneously.

Experimental Set-up
The numerical models were first validated and then calibrated using a published set of experimental data associated with 15 cases. They included 07, 03, and 05 tests with plane slopes, solid submerged structures and porous submerged structures, respectively. The solid submerged structure reference velocity for clear water, discretized to Δ 4 , while the advection terms are differenced wit difference scheme. The computations are carried out for free surface elevati for clear water, , reference velocity for porous medium, and turb simultaneously.

Experimental Set-up
The numerical models were first validated and then calibrated using a publi data associated with 15 cases. They included 07, 03, and 05 tests w submerged structures and porous submerged structures, respectively. The so , reference velocity for porous medium, 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, , reference velocity for porous medium, and turbulent kinetic energy, k simultaneously.

Experimental Set-up
The numerical models were first validated and then calibrated using a published set of experimental data associated with 15 cases. They included 07, 03, and 05 tests with plane slopes, solid and turbulent kinetic energy, k simultaneously. waves, a sponge layer with an adequate width is provided along region 2. Ranasinghe et al. (2009a) artificial energy dissipation term was incorporated when simulating wave transformation over submerged structures to overcome unrealistic flow patterns and numerical instabilities. In addition, considering the bed as permeable with finite 8 Offshore Open Boundary As there is significant wave reflection from submerged breakwaters, it is essential to ab waves, which propagate back to wave incident boundary to prevent a build-up of wave inside the computational domain. This is achieved with the line boundary method proposed et al. (1994). In this method, the neighbouring grid points are divided into two types acco their position with respect to line boundary. Both incident and outgoing waves exist ins boundary (region 1), whereas only outgoing waves exist outside the line boundary (regi shown in Figure 2. To completely dissipate the energy of outgoing waves, a sponge layer adequate width is provided along region 2.  Ranasinghe et al. (2009a) artificial energy dissipation term was incorporated when simulati transformation over submerged structures to overcome unrealistic flow patterns and n instabilities. In addition, considering the bed as permeable with finite porosity, possible wet drying coexisting field (i.e. the wave trough level was allowed to touch the crown of the s was allowed in simulations associated with solid submerged structures following Kenned (2000). However, this was not permitted in simulations associated with porous su structures due to the clear existence of water and porous layers in the model.

Numerical Scheme
The truncated Chen (2006) and Nwogu (1993) Boussinesq-type equations are solved in domain using a third-order Adams-Bashforth predictor step and a fourth-order Adamscorrector step. The computational domain is discretized with grid size, Δ in the -direction equation variables , ℎ, and are defined at grid points in a staggered scheme as show 2. The water depth and surface elevation are defined at grid points …, − 1, , + 1, … velocities are defined half a grid point on either side of elevation grid points. The boundaries of the computational domain correspond to velocity grid points. The first-orde derivatives are differenced to Δ 4 , which automatically eliminates error terms that wou the same form as dispersive terms ]. The second-order spatial deriva

March 2021
Journal of the National Science Foundation of Sri Lanka 49(1)

Experimental set-up
The numerical models were first validated and then calibrated using a published set of experimental data associated with 15 cases. They included 07, 03, and 05 tests with plane slopes, solid submerged structures and porous submerged structures, respectively. The solid submerged structure models were constructed with concrete and the porous were constructed with gravel of mean diameter (d 50 ) equal to 12 mm. Figure 3 illustrates the definition of variables associated with the physical models used and Table 1      Experimental Set-up The numerical models were first validated and then calibrated using a published set o data associated with 15 cases. They included 07, 03, and 05 tests with plane submerged structures and porous submerged structures, respectively. The solid subm models were constructed with concrete and the porous were constructed with g diameter (d50) equal to 12 mm. Figure 3 illustrates the definition of variables asso physical models used and Table 1 lists the wave parameters, physical model d numerical space and time steps associated with all cases considered. A series of numerical simulations were carried out to estimate the best fit for each c the wave breaking detection parameter, / and mixing length, . Here, the wav determined by linear dispersion relation. The wave breaking detection parameter is v equation as it determines the production of turbulent kinetic energy. The / ratio intervals of 0.02 starting from 0.60 to 1.00, while was initially varied as a fracti of 0.25 from 0.25 to 1.50) of the incident wave height making the total number of s for an individual case equal to 126. All the physical and numerical models were test waves only. A series of numerical simulations were carried out to estimate the best fit for each case by varying the wave breaking detection parameter, 17 observed for cases with submerged structures is due to the forced abrupt wave breaking.

Model Calibration
In the present study, using the simulation results, it is attempted to propose appropriate values for wave breaking detection parameter and mixing length to make the model applicable to different wave environments and bottom configurations (i.e. the calibration). The model performance for conditions with plane seabed slopes and submerged structures is considered separately and the wave breaking detection parameter, / through the calibration process is estimated to be 0.90 and 0.80, respectively for these two bottom configurations.
Considering the turbulence generation and associated scale of turbulence, it is hypothesized in this study that the mixing length is dependent on both wave height and wave period (not just limited to wave height). As the wavelength is proportional to the wave period, mixing length, is defined as a fraction of 0 , where 0 is the deepwater wavelength. Table 2 lists this quantity for all 15 cases and the ratio, / 0 is found to be slightly higher for the cases associated with plane seabed slopes in comparison to those associated with submerged structures. The average / 0 ratios are 0.15 and 0.09 for cases with plane seabed slopes and submerged structures, respectively (Table 3).
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 horizontal velocity at a reference water depth, , exceeds a threshold velocity (a fraction of the wave celerity, C). The turbulent length scale, is to be quantified appropriately as part of the study, while 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, and that is obtained by assuming a local balance between production and dissipation of turbulent kinetic energy as follows: . Here, the wave celerity, C is determined by linear dispersion relation. The wave breaking detection parameter is vital in the TKE equation as it determines the production of turbulent kinetic energy. The 17 mixing length, tend to be slightly smaller for cases with submerged structures in comparison to cases with plane seabed slopes. The relatively smaller mixing lengths found for submerged structures are mainly attributed by the depth limitations on the crown and the relatively smaller / observed for cases with submerged structures is due to the forced abrupt wave breaking.

Model Calibration
In the present study, using the simulation results, it is attempted to propose appropriate values for wave breaking detection parameter and mixing length to make the model applicable to different wave environments and bottom configurations (i.e. the calibration). The model performance for conditions with plane seabed slopes and submerged structures is considered separately and the wave breaking detection parameter, / through the calibration process is estimated to be 0.90 and 0.80, respectively for these two bottom configurations.
Considering the turbulence generation and associated scale of turbulence, it is hypothesized in this study that the mixing length is dependent on both wave height and wave period (not just limited to wave height). As the wavelength is proportional to the wave period, mixing length, is defined as a fraction of 0 , where 0 is the deepwater wavelength. Table 2 lists this quantity for all 15 cases and the ratio, / 0 is found to be slightly higher for the cases associated with plane seabed slopes in comparison to those associated with submerged structures. The average / 0 ratios are 0.15 and 0.09 for cases with plane seabed slopes and submerged structures, respectively (Table 3). term for turbulent kinetic energy production by wave breaking, i.e., one equation turbulence closure model (Nwogu (1996)): 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 horizontal velocity at a reference water depth, , exceeds a threshold velocity (a fraction of the wave celerity, C). The turbulent length scale, is to be quantified appropriately as part of the study, while 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, and that is obtained by assuming a local balance between production and dissipation of turbulent kinetic energy as follows: was initially varied as a fraction (in intervals of 0.25 from 0.25 to 1.50) of the incident wave height making the total number of simulation runs for an individual case equal to 126. All the physical and numerical models were tested with regular waves only.

Comparison of measured and simulated quantities
To investigate the performance of the wave breaking induced energy dissipation sub-model, numerical model results are compared with the measured wave height and mean water level distributions across the surf zone. The overall agreement between measured and simulated values at all gauges is quantified firstly by calculating the Willmott index (Wilmott,1984) given in equation 19 for wave height and mean water level separately and then by taking the product between them. 11 numerical model results are compared with the measure distributions across the surf zone. The overall agreement be all gauges is quantified firstly by calculating Willmott inde for wave height and mean water level separately and then b In Equation (19), denotes the measured values, deno average measured values, and is the number of gauges. T 1.0 for a perfect agreement. Figures 4 and 5 depict, respectively, the computed and m mean water level distributions across the surf zone for cas The red lines indicate the profiles associated with the case index (i.e. product between Willmott indices for wave heig the blue lines indicate the profiles associated with the considering the overall agreement of all 07 cases (i.e. usin model parameters and the Willmott indices for cases asso experimental data are listed in Table 2. ... (19) In Equation (19), 11 numerical model results are compared with the meas distributions across the surf zone. The overall agreement all gauges is quantified firstly by calculating Willmott i for wave height and mean water level separately and then

Best Overall Agreement for Individual Cases
In Equation (19), denotes the measured values, de average measured values, and is the number of gauges 1.0 for a perfect agreement. Figures 4 and 5 depict, respectively, the computed an mean water level distributions across the surf zone for The red lines indicate the profiles associated with the c index (i.e. product between Willmott indices for wave h the blue lines indicate the profiles associated with th considering the overall agreement of all 07 cases (i.e. u model parameters and the Willmott indices for cases a experimental data are listed in Table 2. denotes the measured values,

11
To investigate the performance of the wave breaking induced energy dissipation sub-model, numerical model results are compared with the measured wave height and mean water level distributions across the surf zone. The overall agreement between measured and simulated values at all gauges is quantified firstly by calculating Willmott index (Wilmott,1984) given in Equation 19 for wave height and mean water level separately and then by taking the product between them.
In Equation (19), denotes the measured values, denotes the computed values, � denotes the average measured values, and is the number of gauges. The Willmott index will take the value of 1.0 for a perfect agreement. Figures 4 and 5 depict, respectively, the computed and measured wave height distributions and mean water level distributions across the surf zone for cases associated with plane seabed slopes. The red lines indicate the profiles associated with the cases related to highest combined Willmott index (i.e. product between Willmott indices for wave heights and for mean water level), whereas the blue lines indicate the profiles associated with the proposed values for model parameters considering the overall agreement of all 07 cases (i.e. using the calibrated model). The numerical model parameters and the Willmott indices for cases associated with the best agreement with the experimental data are listed in Table 2. denotes the computed values,

11
To investigate the performance of the wave breaking induced energy dissipation sub-model, numerical model results are compared with the measured wave height and mean water level distributions across the surf zone. The overall agreement between measured and simulated values at all gauges is quantified firstly by calculating Willmott index (Wilmott,1984) given in Equation 19 for wave height and mean water level separately and then by taking the product between them.
In Equation (19), denotes the measured values, denotes the computed values, � denotes the average measured values, and is the number of gauges. The Willmott index will take the value of 1.0 for a perfect agreement. Figures 4 and 5 depict, respectively, the computed and measured wave height distributions and mean water level distributions across the surf zone for cases associated with plane seabed slopes. The red lines indicate the profiles associated with the cases related to highest combined Willmott index (i.e. product between Willmott indices for wave heights and for mean water level), whereas the blue lines indicate the profiles associated with the proposed values for model parameters considering the overall agreement of all 07 cases (i.e. using the calibrated model). The numerical model parameters and the Willmott indices for cases associated with the best agreement with the experimental data are listed in Table 2.

Comparison of Measured and Simulated Quantities
To investigate the performance of the wave breaking induc numerical model results are compared with the measured w distributions across the surf zone. The overall agreement betwee all gauges is quantified firstly by calculating Willmott index (W for wave height and mean water level separately and then by taki In Equation (19), denotes the measured values, denotes th average measured values, and is the number of gauges. The W 1.0 for a perfect agreement. Figures 4 and 5 depict, respectively, the computed and measu mean water level distributions across the surf zone for cases as The red lines indicate the profiles associated with the cases rel index (i.e. product between Willmott indices for wave heights a the blue lines indicate the profiles associated with the propo considering the overall agreement of all 07 cases (i.e. using the model parameters and the Willmott indices for cases associated experimental data are listed in Table 2.

Best Overall Agreement for Individual Cases
is the number of gauges. The Willmott index will take the value of 1.0 for a perfect agreement.

Best overall agreement for individual cases
Figures 4 and 5 depict, respectively, the computed and measured wave height distributions and mean water level distributions across the surf zone for cases associated with plane seabed slopes. The red lines indicate the profiles associated with the cases related to highest combined Willmott index (i.e. product between Willmott indices for wave heights and for mean water level), whereas the blue lines indicate the profiles associated with the proposed values for model parameters considering the overall agreement of all 7 cases (i.e. using the calibrated model). The numerical model parameters and the Willmott indices for cases associated with the best agreement with the experimental data are listed in Table 2.
The combined Willmott index for best fit is found to The combined Willmott index for best fit is found to be less than 0.90 only in cases PLS05, PLS06 and PLS07 demonstrating a high-level accuracy of the model across the values used. The wave breaking detection parameter / and mixing length, vary from 0.68 to 1.00, and from 0.75 Hi to 1.50 Hi, respectively across the 07 cases tested with plane seabed slopes. The combined Willmott index for best fit is found to be less than 0.90 only in cases PLS05, PLS06 and PLS07 demonstrating a high-level accuracy of the model across the values used. The wave breaking detection parameter / and mixing length, vary from 0.68 to 1.00, and from 0.75 Hi to 1.50 Hi, respectively across the 07 cases tested with plane seabed slopes.  6 and 7 depict, respectively, the simulated and measured wave height distributions and mean water level distributions across the surf zone for cases associated with submerged structures (both solid and porous). Similar to the cases with plane seabed slopes, the red lines indicate the profiles associated with the cases related to highest combined Willmott index, whereas the blue lines indicate the profiles associated with the proposed values for model parameters considering the overall agreement of all 08 cases. The numerical model parameters and the Willmott indices for cases associated with the best agreement with the experimental data are also listed in Table 2.

, mm
The combined Willmott index for best fit is found to be less than 0.90 only in cases PLS05, PLS06 and PLS07 demonstrating a high-level accuracy of the model across the values used. The wave breaking detection parameter / and mixing length, vary from 0.68 to 1.00, and from 0.75 Hi to 1.50 Hi, respectively across the 07 cases tested with plane seabed slopes. The combined Willmott index for best fit is found to be less than 0.90 only in cases PLS05, PLS06 and PLS07 demonstrating a high-level accuracy of the model across the values used. The wave breaking detection parameter / and mixing length, vary from 0.68 to 1.00, and from 0.75 Hi to 1.50 Hi, respectively across the 07 cases tested with plane seabed slopes. The combined Willmott index for best fit is found to be less than 0.90 only in cases PLS05, PLS06 and PLS07 demonstrating a high-level accuracy of the model across the values used. The wave breaking detection parameter / and mixing length, vary from 0.68 to 1.00, and from 0.75 Hi to 1.50 Hi, respectively across the 07 cases tested with plane seabed slopes.  6 and 7 depict, respectively, the simulated and measured wave height distributions and mean water level distributions across the surf zone for cases associated with submerged structures (both solid and porous). Similar to the cases with plane seabed slopes, the red lines indicate the profiles associated with the cases related to highest combined Willmott index, whereas the blue lines indicate the profiles associated with the proposed values for model parameters considering the overall agreement of all 08 cases. The numerical model parameters and the Willmott indices for cases associated with the best agreement with the experimental data are also listed in Table 2.
The combined Willmott index for best fit is found to be less than 0.90 only in cases PLS05, PLS06 and PLS07 demonstrating a high-level accuracy of the model across the values used. The wave breaking detection parameter / and mixing length, vary from 0.68 to 1.00, and from 0.75 Hi to 1.50 Hi, respectively across the 07 cases tested with plane seabed slopes.  Figures 6 and 7 depict, respectively, the simulated and measured wave height distributions and mean water level distributions across the surf zone for cases associated with submerged structures (both solid and porous). Similar to the cases with plane seabed slopes, the red lines indicate the profiles associated with the cases related to highest combined Willmott index, whereas the blue lines indicate the profiles associated with the proposed values for model parameters considering the overall agreement of all 8 cases. The numerical model parameters and the Willmott indices for cases associated with the best agreement with the experimental data are also listed in Table 2.
The combined Willmott index for best agreement is found to be over 0.90 for all cases with submerged structures demonstrating excellent predictive skills of the numerical model. However, the wave breaking detection parameter and mixing length, vary from 0.68 to 1.00, and from 0.5 H i to 1.00 H i , respectively across the 8 cases tested with submerged structures.
From the results obtained by running the numerical model with different model parameters and different depth configurations, it is evident that both the wave breaking detection parameter and the mixing length, From the results obtained by running the numerical model with different model para different depth configurations, it is evident that both the wave breaking detection param mixing length, tend to be slightly smaller for cases with submerged structures in co cases with plane seabed slopes. The relatively smaller mixing lengths found for structures are mainly attributed by the depth limitations on the crown and the relatively observed for cases with submerged structures is due to the forced abrupt wave breakin

Model Calibration
In the present study, using the simulation results, it is attempted to propose appropriat wave breaking detection parameter and mixing length to make the model applicable tend to be slightly smaller for cases with submerged structures in comparison to cases with plane seabed slopes. The relatively smaller mixing lengths found for submerged structures are mainly attributed by the depth limitations on the crown and the relatively smaller rom the results obtained by running the numerical model with different model parameters and ifferent depth configurations, it is evident that both the wave breaking detection parameter and the ixing length, tend to be slightly smaller for cases with submerged structures in comparison to ases with plane seabed slopes. The relatively smaller mixing lengths found for submerged tructures are mainly attributed by the depth limitations on the crown and the relatively smaller / observed for cases with submerged structures is due to the forced abrupt wave breaking.
.3 Model Calibration n the present study, using the simulation results, it is attempted to propose appropriate values for ave breaking detection parameter and mixing length to make the model applicable to different ave environments and bottom configurations (i.e. the calibration). The model performance for onditions with plane seabed slopes and submerged structures is considered separately and the ave breaking detection parameter, / through the calibration process is estimated to be 0.90 nd 0.80, respectively for these two bottom configurations.
From the results obtai different depth configu mixing length, tend cases with plane seab structures are mainly at observed for cases w

Model Calibra
In the present study, us wave breaking detectio wave environments an conditions with plane wave breaking detectio and 0.80, respectively f Considering the turbule observed for cases with submerged structures is due to the forced abrupt wave breaking.

Model calibration
In the present study, using the simulation results, it is attempted to propose appropriate values for wave breaking detection parameter and mixing length to make

17
observed for cases with submerged structures is due to the forced abrupt wave breaking.

Model Calibration
In the present study, using the simulation results, it is attempted to propose appropriate values for wave breaking detection parameter and mixing length to make the model applicable to different wave environments and bottom configurations (i.e. the calibration). The model performance for conditions with plane seabed slopes and submerged structures is considered separately and the wave breaking detection parameter, / through the calibration process is estimated to be 0.90 and 0.80, respectively for these two bottom configurations.
Considering the turbulence generation and associated scale of turbulence, it is hypothesized in this study that the mixing length is dependent on both wave height and wave period (not just limited to wave height). As the wavelength is proportional to the wave period, mixing length, is defined as a fraction of 0 , where 0 is the deepwater wavelength. Table 2 lists this quantity for all 15 cases and the ratio, / 0 is found to be slightly higher for the cases associated with plane seabed slopes in comparison to those associated with submerged structures. The average / 0 ratios are 0.15 and 0.09 for cases with plane seabed slopes and submerged structures, respectively (Table 3). Considering the turbulence generation and associated scale of turbulence, it is hypothesised in this study that the mixing length is dependent on both wave height and wave period (not just limited to wave height). As the wavelength is proportional to the wave period, mixing length, 17 ons with plane seabed slopes and submerged structures is considered separately and the reaking detection parameter, / through the calibration process is estimated to be 0.90 0, respectively for these two bottom configurations. ering the turbulence generation and associated scale of turbulence, it is hypothesized in this hat the mixing length is dependent on both wave height and wave period (not just limited to eight). As the wavelength is proportional to the wave period, mixing length, is defined as on of 0 , where 0 is the deepwater wavelength. Table 2 lists this quantity for all 15 nd the ratio, / 0 is found to be slightly higher for the cases associated with plane slopes in comparison to those associated with submerged structures. The average / 0 re 0.15 and 0.09 for cases with plane seabed slopes and submerged structures, respectively 3). 3 -Proposed breaking indices, mixing lengths and corresponding Willmott indices ted with all cases.
, mm is defined as a fraction of conditions with plane seabed slopes and sub wave breaking detection parameter, / thro and 0.80, respectively for these two bottom con Considering the turbulence generation and asso study that the mixing length is dependent on b wave height). As the wavelength is proportiona a fraction of 0 , where 0 is the deepwat cases and the ratio, / 0 is found to be seabed slopes in comparison to those associate ratios are 0.15 and 0.09 for cases with plane s (Table 3). conditions with plane seabed slopes and subm wave breaking detection parameter, / throu and 0.80, respectively for these two bottom conf Considering the turbulence generation and asso study that the mixing length is dependent on bo wave height). As the wavelength is proportional a fraction of 0 , where 0 is the deepwate cases and the ratio, / 0 is found to be seabed slopes in comparison to those associated ratios are 0.15 and 0.09 for cases with plane se (Table 3). is the deepwater wavelength. Table 2 lists this quantity  for all 15 cases and the ratio,   17 conditions with plane seabed slopes and submerged struc wave breaking detection parameter, / through the calib and 0.80, respectively for these two bottom configurations.
Considering the turbulence generation and associated scale study that the mixing length is dependent on both wave heig wave height). As the wavelength is proportional to the wave a fraction of 0 , where 0 is the deepwater wavelengt cases and the ratio, / 0 is found to be slightly high seabed slopes in comparison to those associated with subme ratios are 0.15 and 0.09 for cases with plane seabed slopes (Table 3). is found to be slightly higher for the cases associated with plane seabed slopes in comparison to those associated with submerged structures. The average 17 wave environments and bottom configurations (i.e. the calibration). The model performance for conditions with plane seabed slopes and submerged structures is considered separately and the wave breaking detection parameter, / through the calibration process is estimated to be 0.90 and 0.80, respectively for these two bottom configurations.
Considering the turbulence generation and associated scale of turbulence, it is hypothesized in this study that the mixing length is dependent on both wave height and wave period (not just limited to wave height). As the wavelength is proportional to the wave period, mixing length, is defined as a fraction of 0 , where 0 is the deepwater wavelength. Table 2 lists this quantity for all 15 cases and the ratio, / 0 is found to be slightly higher for the cases associated with plane seabed slopes in comparison to those associated with submerged structures. The average / 0 ratios are 0.15 and 0.09 for cases with plane seabed slopes and submerged structures, respectively (Table 3). 14 The combined Willmott index for best fit is found to be less than 0.90 only in cases PLS05, PLS06 and PLS07 demonstrating a high-level accuracy of the model across the values used. The wave breaking detection parameter / and mixing length, vary from 0.68 to 1.00, and from 0.75 Hi to 1.50 Hi, respectively across the 07 cases tested with plane seabed slopes.  Table 2.
The combined Willmott index for best agreement is found to be over 0.90 for all cases with submerged structures demonstrating excellent predictive skills of the numerical model. However,

14
The combined Willmott index for best fit is found to be less than 0.90 only in cases PLS05, PLS06 and PLS07 demonstrating a high-level accuracy of the model across the values used. The wave breaking detection parameter / and mixing length, vary from 0.68 to 1.00, and from 0.75 Hi to 1.50 Hi, respectively across the 07 cases tested with plane seabed slopes.  Table 2.
The combined Willmott index for best agreement is found to be over 0.90 for all cases with submerged structures demonstrating excellent predictive skills of the numerical model. However,

14
The combined Willmott index for best fit is found to be less than 0.90 only in cases PLS05, PLS06 and PLS07 demonstrating a high-level accuracy of the model across the values used. The wave breaking detection parameter / and mixing length, vary from 0.68 to 1.00, and from 0.75 Hi to 1.50 Hi, respectively across the 07 cases tested with plane seabed slopes. The combined Willmott index for best fit is found to be less than 0.90 only in cases PLS05, PLS06 and PLS07 demonstrating a high-level accuracy of the model across the values used. The wave breaking detection parameter / and mixing length, vary from 0.68 to 1.00, and from 0.75 Hi to 1.50 Hi, respectively across the 07 cases tested with plane seabed slopes. Figures 6 and 7 depict, respectively, the simulated and measured wave height distributions and mean water level distributions across the surf zone for cases associated with submerged structures (both solid and porous). Similar to the cases with plane seabed slopes, the red lines indicate the profiles associated with the cases related to highest combined Willmott index, whereas the blue lines indicate the profiles associated with the proposed values for model parameters considering the overall agreement of all 08 cases. The numerical model parameters and the Willmott indices for cases associated with the best agreement with the experimental data are also listed in Table 2.
The combined Willmott index for best agreement is found to be over 0.90 for all cases with submerged structures demonstrating excellent predictive skills of the numerical model. However,

14
The combined Willmott index for best fit is found to be less than 0.90 only in cases PLS05, PLS06 and PLS07 demonstrating a high-level accuracy of the model across the values used. The wave breaking detection parameter / and mixing length, vary from 0.68 to 1.00, and from 0.75 Hi to 1.50 Hi, respectively across the 07 cases tested with plane seabed slopes.  As indicated earlier, the blue lines in Figures 4,5,6,7 indicate the simulated wave height and mean water level distributions after setting the wave breaking detection parameter, As indicated earlier, the blue lines in Figures 4,5,6,7 indicate the simulated wave height a water level distributions after setting the wave breaking detection parameter, / to 0.90 and / 0 ratios to 0.15 and 0.09 for cases with plane seabed slopes and submerged s respectively. As it can be seen from Table 3, the combined Willmott index is found to be 0.90 only in results obtained for cases PLS05, PLS06 and PLS07 with the calibrated mod from the significant deviation observed in case PLS06 for mean water level (Figur simulated wave height and mean water distributions with proposed model parameter val excellent agreement with experimental data. It is noted from the simulations that the mean water levels associated with porous demonstrate a sudden increase over the structure when compared to observations (not p for PSB01 due to its short crest). Although flow through porous structures is allowe simulations, wave celerity calculations were performed using linear dispersion relation tak water depths. Since the porosity effects were not considered, wave celerity might h underestimated in the simulations, which led to a slightly early and intense wave breaking. ratios to 0.15 and 0.09 for cases with plane seabed slopes and submerged structures, respectively. As it can be seen from Table 3, the combined Willmott index is found to be less than 0.90 only in results obtained for cases PLS05, PLS06 and PLS07 with the calibrated model. Apart from the significant deviation observed in case PLS06 for mean water level (Figure 5), the simulated wave height and mean water distributions with proposed model parameter values show excellent agreement with experimental data.
It is noted from the simulations that the mean water levels associated with porous structures demonstrate a sudden increase over the structure when compared to observations (not prominent for PSB01 due to its short crest). Although flow through porous structures is allowed in the simulations, wave celerity calculations were performed using linear dispersion relation taking clear water depths. Since the porosity effects were not considered, wave celerity might have been Journal of the National Science Foundation of Sri Lanka 49 (1) March 2021 underestimated in the simulations, which led to a slightly early and intense wave breaking.

CONCLUSIONS
Two weakly nonlinear Boussinesq-type wave models were developed, validated and calibrated using a published set of data in 1DH. During the calibration process, wave breaking detection parameter, 19

Conclusions
Two weakly nonlinear Boussinesq-type wave models are developed, validated and calibrated using a published set of data in 1DH. During the calibration process, the wave breaking detection parameter, / and mixing length, were varied obtaining a total of 126 simulation runs fo each experimental case considered (there were 15 cases). Using the simulation results in this study the mixing length is correlated to both wave height and wavelength considering its direct relevance to those wave parameters and is defined as a fraction of 0 . The simulation runs indicate a parametric value of 0.90 and 0.80 for / , and a non-dimensional mixing length ( / 0 ) o 0.15 and 0.09, for modelling of wave transformation over plane seabeds and submerged structures respectively. The results obtained from 1DH numerical experiments show promising results for it capabilities in simulating wave evolution over different depth configurations in 2DH. and mixing length,

19
lusions eakly nonlinear Boussinesq-type wave models are developed, validated and calibrated using shed set of data in 1DH. During the calibration process, the wave breaking detection ter, / and mixing length, were varied obtaining a total of 126 simulation runs for perimental case considered (there were 15 cases). Using the simulation results in this study, ing length is correlated to both wave height and wavelength considering its direct relevance e wave parameters and is defined as a fraction of 0 . The simulation runs indicate a tric value of 0.90 and 0.80 for / , and a non-dimensional mixing length ( / 0 ) of d 0.09, for modelling of wave transformation over plane seabeds and submerged structures, ively. The results obtained from 1DH numerical experiments show promising results for its ities in simulating wave evolution over different depth configurations in 2DH.
were varied obtaining a total of 126 simulation runs for each experimental case considered (there were 15 cases). Using the simulation results in this study, the mixing length is correlated to both wave height and wavelength considering its direct relevance to those wave parameters and is defined as a fraction of 19 4. Conclusions Two weakly nonlinear Boussinesq-type wave models are developed, validated and calibrated using a published set of data in 1DH. During the calibration process, the wave breaking detection parameter, / and mixing length, were varied obtaining a total of 126 simulation runs for each experimental case considered (there were 15 cases). Using the simulation results in this study, the mixing length is correlated to both wave height and wavelength considering its direct relevance to those wave parameters and is defined as a fraction of 0 . The simulation runs indicate a parametric value of 0.90 and 0.80 for / , and a non-dimensional mixing length ( / 0 ) of 0.15 and 0.09, for modelling of wave transformation over plane seabeds and submerged structures, respectively. The results obtained from 1DH numerical experiments show promising results for its capabilities in simulating wave evolution over different depth configurations in 2DH.
The simulation runs indicate a parametric value of 0.90 and 0.80 for 19 clusions eakly nonlinear Boussinesq-type wave models are developed, validated and calibrated using lished set of data in 1DH. During the calibration process, the wave breaking detection eter, / and mixing length, were varied obtaining a total of 126 simulation runs for xperimental case considered (there were 15 cases). Using the simulation results in this study, ixing length is correlated to both wave height and wavelength considering its direct relevance se wave parameters and is defined as a fraction of 0 . The simulation runs indicate a etric value of 0.90 and 0.80 for / , and a non-dimensional mixing length ( / 0 ) of nd 0.09, for modelling of wave transformation over plane seabeds and submerged structures, tively. The results obtained from 1DH numerical experiments show promising results for its ilities in simulating wave evolution over different depth configurations in 2DH.
, and a non-dimensional mixing length loped, validated and calibrated using ocess, the wave breaking detection g a total of 126 simulation runs for g the simulation results in this study, ength considering its direct relevance 0 . The simulation runs indicate a sional mixing length ( / 0 ) of e seabeds and submerged structures, ments show promising results for its onfigurations in 2DH.
of 0.15 and 0.09, for modelling of wave transformation over plane seabeds and submerged structures, respectively. The results obtained from 1DH numerical experiments show promising results for its capabilities in simulating wave evolution over different depth configurations in 2DH.