Energy transfer due to shoaling and decomposition of breaking and non-breaking waves over a submerged bar

ABSTRACT Wave propagation over a submerged bar is simulated using the open source CFD model REEF3D with various incident wave heights to study shoaling, wave breaking features and the process of wave decomposition into higher harmonics for relatively long waves of kd=0.52. The computed free surface elevations are compared with experimental data and good agreement is obtained for both non-breaking and spilling breaking waves for both the wave phase and free surface elevation, which has been difficult to obtain in current literature. The differences in the mode of wave shoaling over the weather side slope and the wave decomposition over the lee side slope of the submerged bar are discussed. The evolution of spilling breakers and plunging breakers over the bar crest is also studied. It is found that the free surface elevation continuously increases due to shoaling in the case of non-breaking waves, whereas breaking waves propagate with much lower free surface elevations after breaking over the bar crest. The power spectra of the free surface elevations at various locations indicate that the wave energy in the fundamental frequency is reduced by 76 for the non-breaking wave with kA=0.015 and by about 90 in other cases with higher incident wave heights with kA=0.023−0.034 due to energy dissipation and energy transfer to higher harmonic components as the wave propagates over the submerged bar.


Introduction
Wave propagation in shallow waters is influenced by sea bottom topography and wave transformation processes such as diffraction, shoaling and wave breaking. Wave shoaling refers to the phenomenon where the incident wave height is changed as the deep water wave propagates to water depths less than half the wavelength. Shoaling results in asymmetry in the wave profile with sharper crests and shallower troughs, creating an imbalance in the local wave energy distribution and leading to wave deformation (Adeyemo, 1969). The wave crest heights reach a limiting value, beyond which the wave breaks to balance the local increase in the wave energy. The additional challenge in wave propagation over a submerged obstacle is the wave decomposition process which occurs behind the obstacle, in the region of increasing water depth, leading to the evolution of higher-order harmonics and rapidly varying waveforms. These processes can only be represented in a numerical model that accounts for nonlinearity and has good dispersion characteristics (Beji & Battjes, 1994).
The accurate evaluation of the wave kinematics in the near-shore area is important due to their impact on CONTACT Arun Kamath arun.kamath@ntnu.no hydrodynamic properties such as wave forces, wave runup and sediment transport. The mode of wave breaking is generally classified using the surf similarity parameter, ξ = tan α/ √ H/L 0 , where α is the angle of the slope, H is the incident wave height over the toe of the slope and L 0 = (g/2π)T 2 is the deep water wavelength, where T is the wave period. Battjes (1974) presented the relationships between ξ and various flow parameters, and also the classification of breaker types on emergent plane sloping beaches. Gourlay (1994) carried out experiments on waves breaking on a submerged reef and Blenkinsopp and Chaplin (2008) on a submerged slope and found that the classification presented by Battjes (1974) for emergent sloping beaches is not directly applicable for submerged structures. Wave propagation over submerged structures has been studied through experimental investigations on a submerged bar (Beji & Battjes, 1993), a rectangular obstacle (Chang, Hsu, & Liu, 2001) and processes such as wave decomposition and vortex generation have been identified. Numerical modeling of wave propagation over a submerged obstacle has been carried out using Boussinesq equations (Beji & Battjes, 1994;Bosboom, Klopman, Roelvink, & Battjes, 1996;Brocchini, Drago, & Iovenitti, 1992) and shallow water equations (Kobayashi, Otta, & Roy, 1987) with good results for the wave shoaling process. Lemos (1992) was the first to present simulations of breaking waves using the Reynolds-Averaged Navier-Stokes (RANS) equations. Lin and Li (1998) and Zhao, Armfield, and Tanimoto (2004) employed singlephase Computational Fluid Dynamics (CFD) models to simulate breaking waves, which could not account for the air-water interaction responsible for the complex free surface deformations (Christensen, 1998).
The knowledge of wave transformation and transmission across submerged structures finds its application in coastal protection measures such as submerged breakwaters, ecological conservation and recreational measures such as artificial reefs. The wave decomposition process modifies the waves transmitted over the submerged structure and this can be exploited usefully in a combined submerged bar-floating breakwater coastal protection measure. For the design of recreational artificial reefs and bars, it is essential to have a better idea regarding the breaking wave characteristics on the crest of the bar to provide sufficient breaker heights for surfing. It has been presented in previous studies on wave breaking that the wave breaking characteristics vary significantly under different breaking conditions (Battjes, 1974;Blenkinsopp & Chaplin, 2008;Gourlay, 1994). In addition, the many existing numerical and theoretical models for wave transformation over submerged breakwaters are based on the potential flow assumption, which cannot describe the rotational flow that occurs during the breaking process (Brocchini, 2013). Numerical modeling of shallow water flows (Lai & Khan, 2012;Muttil & Chau, 2007), in estuaries (Chau & Jiang, 2001Wu, & Chau, 2006) and over hydraulic structures (Haun, Olsen, & Feurich, 2011) have been numerically investigated by several authors. CFD modeling solves the fluid flow problem by solving the Navier-Stokes equations, accounting for most of the fluid physics with few assumptions. This method has been previously applied to the simulation of breaking waves over a slope by Hieu, Katsutoshi and Ca (2004) and Jacobsen, Fuhrman and Fredsøe (2012) using a Volume of Fluid (VoF)-based interface capturing method and Alagan Chella, Bihs, Myrhaug, and Muskulus (2015) using the level-set method to obtain the interface. Alagan  obtained good agreement with experimental data, with a sharp representation of the breaking wave and the formation of air pockets, due to the higherorder discretization schemes used in the model along with the level-set method, compared with the lowerorder schemes used in previous studies. Numerical modeling with a two-phase CFD model can account for the wave breaking process with few assumptions. Along with higher-order discretization schemes and sharp interface capturing, it can account for the complex free surface process involved during wave transformation including breaking and decomposition.
Several studies have been carried out using Large Eddy Simulation (LES) representation of the turbulence in the Navier-Stokes equations such as those by Christensen & Deigaard (2001), Hi et al., (2004, Zhao et al., (2004), Watanabe, Saeki, and Hosking (2005), Christensen (2006), Lubin, Vincent, Abadie, and Caltagirone (2006) and Mo, Jensen and Liu (2013) to study the wave breaking process. These studies, however, are restricted to studies of the vortex structures under breaking waves, related aspects of turbulence and wave breaking on slopes, but do not represent the wave transformation and decomposition process involved in wave propagation over a submerged bar. While the LES and Direct Numerical Simulation (DNS) approach to CFD modeling of breaking waves can provide a thorough understanding of the micro-scale turbulent features, they are rather expensive and superfluous for engineering applications. The energy transfer and wave decomposition aspects are important from an engineering point of view to determine the efficiency of submerged structures to dissipate wave energy, which can be well resolved using Reynolds-Averaged Navier-Stokes (RANS) models. With a better understanding of this phenomenon, submerged structures can be used as effective coastal protection measures without causing a large visual impact to the coastline, unlike that from an emerged rubble mound breakwater.
In the current study, the open source CFD model REEF3D (Bihs, Kamath, Alagan Chella, Aggarwal, & Arntsen, 2016) is used to simulate wave propagation over a submerged bar. The numerical results are compared with the experimental data from Beji and Battjes (1993). Several previous studies regarding this have numerically calculated the wave propagation only for the non-breaking wave cases (Morgan et al., 2010;Roeber, Cheung, & Kobayashi, 2010;Stelling, & Zijlema, 2003). The breaking wave case was modeled by Tissier, Bonneton, Marche, Chazel, and Lannes (2012), but they reported deviations from the experimental observations from the point of wave breaking. Thus, numerical models accounting for both breaking and non-breaking waves over a submerged bar (Beji & Battjes, 1993) and having good agreement with experimental data for both the free surface elevation and the wave phase have not been presented in the current literature. This is especially true for longer waves with T = 2.5 s, where the wave decomposition process is seen to be much stronger in experiments compared with shorter waves with T = 1 s (Beji & Battjes, 1994). An initial study for only non-breaking wave shoaling on a submerged bar was presented, together with a comparison with experimental data, by Kamath, Alagan Chella, Bihs, and Arntsen (2015).
In this paper, the study is significantly extended to cover the evolution of spilling and plunging breakers on the bar crest, with a comparison of the free surface elevation in the spilling case. In addition, the shoaling process for the different incident waves is examined through a comparison of the relative wave crest elevations and an evaluation of the maximum wave crest steepness and the relative phase differences between the primary wave crests of the transformed waves in the different cases. The decomposition process with transfer of wave energy to higher harmonics is also examined using the power spectral density computed from the free surface elevations, and the redistribution of the wave energy amongst the harmonics is discussed. The effect of wave breaking on the wave propagation and decomposition process is also discussed. The study also addresses the effect of the lee side slope of the submerged bar on the wave decomposition process. This is essential, as the current literature has dealt with the details of the effect of the weather side slope on the wave shoaling and breaking process, but not the details of the wave decomposition on the lee side. This is important for the construction of innovative coastal and harbour protection structures such as a combination of a submerged bar with a floating breakwater. The floating breakwater is efficient at absorbing wavelengths in the order of its width, but not longer waves. The addition of a submerged bar in front of a floating breakwater can be used to dissipate some wave energy and also decompose the waves into shorter waves that can be effectively absorbed by the floating breakwater.

Governing equations
The numerical model uses the Reynolds-Averaged Navier-Stokes (RANS) equations along with the continuity equation to evaluate the fluid flow: where u i is the ensemble averaged velocity, ρ is the density of the fluid, p is the modified pressure, ν is the kinematic viscosity, ν t = k/ω is the eddy viscosity, k is the turbulent kinetic energy, ω is the specific turbulent dissipation rate, t is time and g is the acceleration due to gravity. The equations are presented in the compact tensor notation and indices i and j are free indices. Surface tension is included in the source term S in Equation (2), where where σ is the surface tension coefficient, κ is the curvature of the interface, the surface tension is activated only around the interface using the Dirac delta function δ(φ) (Peng, Merriman, Osher, Zhao, & Kang, 1999) and φ is the level-set function used to obtain the free surface in the model. The projection method (Chorin, 1968) is used for pressure treatment and the resulting Poisson pressure equation is solved using a preconditioned BiCGStab solver (van der Vorst, 1992). Turbulence modeling is carried out using the two-equation k-ω model proposed by Wilcox (1994). Wave propagation is characterized by large gradients in the velocities resulting in a highly strained flow. Since the production of turbulence in the k-ω model depends on the gradients in the velocity field, this results in unphysical overproduction of turbulence in a numerical wave tank. A stress limiter in the definition of eddy viscosity using an assumption of Bradshaw, Ferriss, and Atwell (1967) as shown by Durbin (2009) is implemented to avoid this. The free surface is a natural boundary for the turbulent eddies, which is not accounted for in the k-ω model, resulting in an overproduction of turbulence at the free surface in a two-phase CFD model due to the large strain caused by the large difference in the densities of air and water. Free surface turbulence damping using a limiter around the interface as shown by Naot and Rod (1982) is carried out to avoid the overproduction of turbulence at the interface. The limiter is activated only around the interface using the Dirac delta function.

Discretization schemes
The fifth-order conservative finite difference Weighted Essentially Non-Oscillatory (WENO) scheme (Jiang & Shu, 1996) is used for the discretization of the convective terms in the RANS equations and the level-set function, the turbulent kinetic energy and the specific turbulence dissipation rate are discretized using the Hamilton-Jacobi formulation of the WENO scheme (Jiang & Peng, 2000). Time advancement is carried out using a four-step scheme proposed by Choi and Moin (1994) with implicit treatment of convective and viscous terms. An adaptive time stepping approach is used to satisfy the Courant-Frederich-Lewy (CFL) condition for numerical stability. The numerical model uses a uniform Cartesian grid for spatial discretization facilitating an easy implementation of higher-order schemes. The staggered grid approach is used with pressure at the cell centers and velocities at the cell faces, providing a tight coupling between the pressure and the velocity. A local directional ghost cell immersed boundary method (Berthelsen & Faltinsen, 2008) extended to three dimensions is employed to handle complex geometries. The numerical model is completely parallelized using the MPI library and can be executed on high performance computing systems.

Free surface
The free surface determines the boundary between the two fluids in the numerical model. The distinguishing characteristics of the two fluids are density and viscosity. The free surface in the numerical wave tank is captured using the level-set method (Osher & Sethian, 1988) and the respective material properties are applied to the two phases. Using the level-set method, the interface is represented by the zero level set of the signed distance level-set function φ. The level-set function provides the least distance of each point in the domain from the interface. The different fluids are distinguished by the sign of the level-set function as shown in Equation (4): The definition of the level-set function makes it smooth across the interface and provides a sharp representation of the interface. The level-set function is convected by the velocity field in the numerical wave tank. The signed distance property is lost on convection and is restored by re-initializing the level-set function after every iteration with the partial differential equation re-initialization procedure due to Peng et al. (1999).

Numerical wave tank
The numerical wave tank uses the relaxation method (Larsen & Dancy, 1983) for wave generation and absorption. In this method, relaxation functions are used to moderate the computational values to the expected values from wave theory to generate and absorb waves. This requires certain zones of the wave tank to be reserved as relaxation zones for wave generation and absorption. The numerical model uses the relaxation functions proposed by Jacobsen et al. (2012) presented in Equation (5): where (x) is the relaxation function and x is the coordinate along the x-axis scaled to the length of the relaxation zone. The relaxation functions prescribe the required values for free surface elevation and velocity from wave theory to the wave tank using Equation (6): The relaxation function also absorbs reflections from the objects placed in the numerical wave tank, so that it does not affect wave generation and simulates a wavemaker with active absorption. At the numerical beach, the computational values from the wave tank are reduced to zero so as to absorb the wave energy smoothly without spurious reflections from the beach. A no-slip boundary condition is applied on the bottom wall and on the surface of the objects in the tank and symmetry boundary conditions on the top of the numerical wave tank. The boundary conditons are enforced through the ghost cell immersed boundary method.

Results
A grid refinement study is carried out first to select the grid size to be used for the simulations in the study. Then, wave propagation over a submerged bar is simulated for different incident wave heights and the numerical results are compared with experimental data. The wave transformation over the submerged bar is studied using the data obtained from wave gages at different locations along the length of the bar. The evolution of spilling and plunging breakers on the bar crest in the simulation is also presented. The shoaling process for the different incident waves is examined through the variation of the relative wave crest elevations. The decomposition process with transfer of wave energy to higher harmonics is examined by calculating the power spectral densities of the computed free surface elevations at the different locations in the wave tank.

Grid refinement study
Accurate wave generation and propagation in the numerical wave tank is verified by carrying out a grid refinement study.  from a grid size of dx = 0.02 m onwards. Owing to the high-order discretization schemes used in the model and the relatively low wave steepness in the study, there is no significant difference in the wave heights obtained at the different grid sizes. But, in order to capture the evolution of wave shoaling and breaking in this study, a grid size of dx = 0.005 m is used for the simulations.
In addition to the grid refinement study for wave generation and propagation in the numerical wave tank without any obstacles, a grid refinement study for the plunging breaking wave obtained for H 4 = 0.052 m is carried out with simulations at grid sizes dx = 0.005, 0.01, 0.02 and 0.04 m. This is done to verify that the grid resolution is sufficient for the representation of the complex case of wave breaking. From the free surface profiles in Figure 2, it can be concluded that the breaker location has converged to x = 17.2 m from dx = 0.01 m onwards, but the vertical profile of the breaking wave crest is best represented by dx = 0.005 m. This confirms that the choice of dx = 0.005 m as the grid size for the simulations is justified. The very fine grid required to represent the wave breaking in this study arises from the fact that the incident waves are of low steepness and they undergo large and rapid changes in their wave steepness during propagation over the bar. This follows the conclusions by Alagan  that incident waves with lower steepnesses undergo larger deformations than waves with larger incident steepnesses. The computational time required for this fine grid simulation was about 40 h on 128 processors for a simulation length of 60 s.

Numerical wave tank setup
The simulations of wave propagation over a submerged bar are carried out based on the experimental studies of Beji and Battjes (1993). The submerged bar has a weather side slope of 1 : 20, a lee side slope of 1 : 10 and a crest height of 0. Wave gages (WG1-WG8) are placed at various locations along the bar to evaluate the wave propagation over the bar as shown in Figure 3. A two-dimensional numerical wave tank 38 m long and 0.8 m high with a grid size of 0.005 m is used, resulting in a total of 1.216 million cells. A CFL number of 0.1 is used. A wave generation relaxation zone of length 5 m and a numerical beach of length 9.5 m are used at the beginning and end of the wave tank, respectively, to ensure good wave generation and absorption. The x-coordinate in the wave tank begins at the end of the wave generation relaxation zone and the same distances as in the experiments by Beji and Battjes (1993)

Non-breaking wave propagation over a submerged bar
A simulation is carried out with second-order Stokes waves of wave height H 1 = 0.022 m, wave period T = 2.5 s and wavelength L = 4.74 m. The free surface elevations are computed at several locations along the submerged bar and are compared with the measured experimental data in Figure 4. Good agreement is seen in both the phase and amplitude of the transformed waves.
As the waves propagate along the reducing water depth along the upward slope of the bar, the wave profile is seen to be slightly deformed with the development of a saw-toothed profile at x = 11.0 m in Figure Figure 4, which becomes prominent at x = 12.0 m in Figure 4(c). As a result of wave shoaling, high and sharp wave crests are formed over the bar crest x = 13.0 m in Figure 4(d).
The decomposition of the wave with the development of higher harmonic components is also observed from x = 14.0 m onwards (Figure 4(e)), as the wave propagates over the end of the bar crest. As the wave propagates along the lee side slope of the bar, the water depth increases and a process opposite to wave shoaling takes place as described by Beji and Battjes (1993). The free surface elevation begins to reduce compared with the elevations on the upward slope and the crest. The wave decomposition results in the formation of secondary and tertiary waves after the bar crest as seen in Figures 4(f), 4(g) and 4(h).

Breaking wave propagation over a submerged bar
The incident wave height is further increased to H 3 = 0.042 m to simulate spilling breakers and the computed free surface elevations are compared with experimental  The evolution of the wave profile in the region of wave breaking in the simulation is presented in Figure 6 to obtain further insight into the breaking process in this case. The shoaling of the wave due to the reducing water depth leads to a sharp wave crest on the bar crest as seen in Figure 6(a). The bar crest acts as a flat bottom with very low water depth and the wave propagates over the crest without much change to its amplitude, but with reduced wave celerity. The reduction in wave celerity combined with an increase in wave crest elevation due to shoaling leads to a local imbalance in the wave energy with wave crest particle velocities higher than the wave celerity. This increases the asymmetry of the wave and the appearance of a steep wave crest. The steep wave crest then stretches away from the main wave crest in Figure 6(b). Due to lack of further excess energy, the wave crest then begins to spill forward onto the main wave crest in Figures 6(c) and 6(d), resulting in a small-scale spilling breaker. The velocity contours in the figures demonstrate the large increase in the horizontal water particle velocity in the overturning crest compared with the rest of the free surface, signifying the complex hydrodynamics involved in the spilling breaking wave. The total duration from the near vertical wave crest profile until the wave crest rejoins the preceding trough is only 0.1T, signifying the rapid and small-scale nature of the spilling breaker in this case. Similarly the evolution of the plunging breaking wave over the bar crest for H 4 = 0.052 m is shown in Figures  7(a)-7(e).

Wave transformation process
The variation in the relative wave crest elevations computed at the different wave gages is studied to gain a comparative perspective of the wave transformation process for both non-breaking and breaking waves. The incident waves at x = 6.0 m in Figure 8  small horizontal asymmetry in the wave profile with shallower troughs and sharper crests, which is characteristic of second-order Stokes waves. The breaking and nonbreaking waves show certain differences in the transformation properties. In the case of the non-breaking waves with H 1 = 0.022 m and H 2 = 0.035 m, shoaling leads to saw-toothed asymmetry in the wave profile. The higher incident wave H 2 undergoes a greater increase in the relative crest elevation and attains a sharper sawtoothed asymmetry in Figure 8(b) at x = 11.0 m. As the wave reaches the crest of the bar at x = 12.0 m, the relative crest elevation is higher for H 2 compared with H 1 in Figure 8  & Myrhaug, 1978) = η /L , where η is the wave crest height and L the distance from the wave crest to the wave zero-crossing location, can be used to quantify the crest steepness. In the case of H 2 = 0.035 m, the maximum wave crest steepness max = 0.2008 is calculated for WG4 at x = 13.0 m. The maximum wave crest steepness max = 0.1275 for H 1 = 0.022 m is obtained at x = 14.0 m at WG5. Also, the higher incident wave (H 2 ) moves faster than the lower incident wave (H 1 ). This follows from shallow water wave propagation, where a higher wave propagates faster for a given wave period and water depth. The higher incident wave attains the highest crest elevation during its propagation over the upward slope and thus propagates faster over the shallow water depth over the crest.
The submerged bar crest ends at x = 14.0 m and the initiation of wave decomposition is seen in Figure 8(e), with the appearance of secondary crests. As the wave propagates further, the water depth increases over the downward slope of the submerged bar. This change in the water depth begins a process of de-shoaling (Beji & Battjes, 1993), where the waves reduce in amplitude as they propagate over gradually increasing water depths. Well-developed secondary wave crests are seen at x = 15.0 m in Figure 8(f). It is also observed that the reduction in the relative crest elevation is lower for the lower incident wave. The higher non-breaking wave H 2 , which had the highest crest elevation at x = 13.0 m has a lower primary crest elevation and a higher secondary crest elevation compared with H 1 , indicating that H 2 transfers a larger amount of wave energy to a higher frequency. In Figures 8(g) and 8(h), the formation of a secondary and a tertiary wave crest is seen for both the non-breaking waves. The lower wave (H 1 ) continues to maintain a higher primary relative crest elevation throughout the wave decomposition process, whereas the higher wave has slightly higher secondary and tertiary relative crest elevations.
In  Figures 8(d)-8(h). This is justified by the fact that H 4 evolves into a plunging breaking wave and dissipates a larger part of its energy in the process compared with the spilling breaking wave formed by H 3 . It is also noticed that the surf similarity number is ξ = 1.0623 for H 3 and ξ = 0.9547 for H 4 . According to the classification of Battjes (1974), these correspond to plunging wave breaking on an emergent plane slope, but the results in this study show a spilling breaker for H 3 and a plunging breaker for H 4 . This indicates that the wave breaking on a bar crest has different breaker characteristics and the original classification for wave breaking on emergent plane slopes cannot be directly applied to wave breaking over a submerged bar.
In order to understand the wave transformation process further for the different cases simulated, the phase difference between the different waves during their propagation over the bar is analyzed at the various gage locations. The relative phase difference δθ is calculated as where δx is the horizontal distance between the two primary wave crests and L is the wavelength. The relative phase difference δθ between the primary wave crests for H 2 = 0.035 m, H 3 = 0.042 and H 4 = 0.052 m with respect to H 1 = 0.022 m is presented in Figure 9. It can be concluded that the higher waves As the waves propagate further to x = 17.0 m (WG8), the phase differences return to the values obtained at the end of the bar crest at x = 14.0 m (WG5). Thus, the wave transformation process for all the incident waves is similar up to the region of wave breaking, with a higher incident wave attaining a higher relative crest elevation. After the region of wave breaking, the transformation of the breaking waves depends on the type of wave breaking, whereas the non-breaking waves continue with the trend seen on the weather side slope. In the region of increasing water depth after the bar crest, the lower nonbreaking wave maintains a higher primary relative crest elevation compared with the secondary relative crest elevation. The breaking waves show similar relative crest amplitudes. The higher incident waves are also seen to propagate faster and increase in celerity over the bar up to the end of the bar crest.

Wave decomposition process
The wave decomposition process is examined by calculating the power spectral densities from the free surface elevations at the different locations for the different cases presented in the study. The location of the peaks in the spectra are used to identify the harmonics to which the wave energy is transferred during wave decomposition. The free surface elevations are provided to the Fast Fourier Transform (FFT) algorithm with a sampling interval of t s = 0.01 s over a simulation period of 60 s resulting in 6000 FFT points. Zero-padding is applied to obtain total input points to the closest power of 2, that is, 8192 points. The Nyquist frequency is 50 Hz. The normalized power spectra at the different wave gage locations along the bar are presented in Figure 10. The power spectra for all cases are normalized with the spectral amplitude at the primary wave frequency, f 0 = 0.4 Hz, S max in Figure 10(a). The process of shoaling results in an increase of the energy content at the primary frequency and the first harmonic of the non-breaking waves. The lowest incident wave, H 1 , adds 0.25S max to its fundamental frequency f 0 , whereas H 2 and H 3 gain 0.12S max and 0.16S max , respectively, at the first harmonic at x = 11.0 m. On the other hand, the spectral power density for the highest incident wave H 4 is reduce by 0.05S max at f 0 and increased by 0.14S max at the first harmonic f 1 . At x = 12.0 m, the waves reach the bar crest, and significant spectral densities are obtained up to the fourth harmonic f 4 with 0.024S max for H 4 . Wave breaking occurs between x = 12.0 m and x = 13.0 m for H 3 and H 4 . This corresponds with a reduction of the spectral power density to 0.368S max for H 3 and 0.272S max for H 4 at f 0 at x = 13.0 m. In the case of the non-breaking waves, H 1 retains 0.82S max at f 0 and transfers 0.20S max and 0.11S max to f 1 and f 2 , respectively. As the waves propagate across the bar crest and in the region of increasing water depth at x = 15.0 m, the major portion of the energy is distributed between f 0 and f 2 for H 1 . For H 2 , H 3 and H 4 the major portion of the energy is distributed amongst f 0 , f 2 and f 3 at x = 15.0 m in Figure 10  The variation of the spectral power density in the first four harmonics over the submerged bar for all the four waves is presented in Figure 11. The power spectral density at f 0 is reduced significantly at x = 17.0 m to 0.24S max , 0.11S max , 0.09S max and 0.09S max for H 1 , H 2 , H 3 and H 4 , respectively, in Figure 11(a). From Figure 11(b), it is seen that the first harmonic f 1 initially gains energy in all cases, but loses this energy gradually in all cases except H 1 .
The second harmonic f 2 gradually gains energy as the waves propagate over the bar, with a maximum of 0.42S max at x = 16.0 m for H 1 in Figure 11(c). The maximum spectral power in the third harmonic f 3 appears between x = 14.0 m and x = 15.0 m in Figure 11(d). The following distinct pattern emerges regarding the energy transfer between the different harmonics. The fundamental frequency gradually loses most of its energy as the wave propagates over the bar. The first harmonic gains energy initially on the weather side slope, but loses this energy gradually. The second harmonic gains energy steadily and holds most of the wave energy towards the end of the bar. The third harmonic contains significant amounts of energy in the intermediate stages between x = 13.0 m and x = 15.0 m. Finally, the variation of the total energy in the first four harmonics over the length of the bar for all cases is presented in Figure 12. It is clear that all the waves lose significant amounts of energy except for H 1 , with total spectral power densities of 0.89S max , 0.49S max , 0.31S max and 0.30S max for H 1 , H 2 , H 3 and H 4 , respectively, at x = 17.0 m, losing 1.25S max , 0.53S max , 0.73S max and 0.85S max during propagation over the bar.

Effect of lee side slope on wave decomposition
After a thorough investigation of the wave transformation and decomposition process over a submerged bar, the effect of the lee side slope on the wave decomposition process is studied for a non-breaking wave with H 2 = 0.035 m, a spilling breaking wave with H 3 = 0.042 m and a plunging breaking wave with H 4 = 0.052 m. Additional simulations are carried out using leeside slopes of 1 : 5 and 1 : 20 and compared with the results obtained for the lee side slope of 1 : 10 in the previous section. The free surface elevations are measured at relative distances x r = x g − x t /L, where x g is the gage location, x t is the location of the leeside toe and L is the wavelength. The free surface elevations at x r = −0.10, 0.0, 0.21 and 0.84 are measured and the spectral power density is calculated. The normalized power density spectra at various locations over the submerged bar for H 2 = 0.035 m are presented in Figure 13. It is seen that most of the wave energy is present in the second harmonic f /f 0 = 3 for all the three lee side slopes. The wave energy in the fundamental frequency f /f 0 = 1 is similar for all the three slopes over the region after the submerged bar. In the case of the steeper lee side slope of m = 1 : 5, significant energy is  contained in the third harmonic f /f 0 = 4 at x r = −0.10 and this energy is dissipated as the wave propagates to x r = 0.84. On the other hand, for the milder slope of m = 1 : 20, the wave energy is mostly concentrated in the second harmonic f /f 0 = 3 and is reduced by 50% as the wave propagates from x r = −0.10 to 0.84. The waves propagating over the mildest slope with m = 1 : 20 generally have the lowest energy in the three slopes studied. Thus, for non-breaking waves, a steeper slope leads to energy transfer to the second and third harmonics, and a milder slope results in higher energy in the second harmonic.
In the case of the spilling breaking waves with incident wave height H 3 = 0.042 m, the normalized power spectral density is calculated using the free surface elevations at x r = −0.10, 0.0, 0.21 and 0.84 and the results are presented in Figure 14. Owing to the wave breaking process, the wave energy at x r = −0.10 is lower in this case in Figure 14(a) compared with the non-breaking waves, but the pattern of the wave energy distribution amongst the harmonics is similar. The wave energy for both slopes m = 1 : 5 and 1 : 20 is lower than that for m = 1 : 10 at x r = −0.10. Over the lee side toe of the bar at x r = 0.0, the energy distribution for m = 1 : 5 and m = 1 : 10 is similar up to the second harmonic, whereas the steeper slope allows for more energy in the third harmonic. In the case of the milder slope m = 1 : 20, the wave energy is the lowest amongst the three over all the harmonics except the first harmonic. On propagating further, the wave energy in the fundamental frequency is retained at similar values for all the three slopes, whereas the wave energy in the higher harmonics is reduced significantly, see Figures 14(c) and 14(d). Thus, for spilling breaking waves, the lee side slope does not have a strong influence on the wave transformation and decomposition process after wave breaking. The maximum energy is contained in the second harmonic for all the three slopes studied.
The normalized spectral power density at different locations over a submerged bar for plunging breaking waves propagating over different lee side slopes is presented in Figure 15. Similar to the spilling breaking waves, the lee side slope of m = 1 : 10 has the highest energy in the third harmonic f /f 0 = 3. The waves propagating over the steepest slope of m = 1 : 5 have the highest energy in the fundamental frequency f /f 0 = 1. Thus, a steeper lee side slope results in higher wave energy being retained in the fundamental and second harmonic while a milder slope results in lower energy in the fundamental frequency. The main difference between the case with spilling breaking and plunging breaking is the higher total energy loss in the case of plunging braking waves.
From the results for the wave decomposition process for non-breaking, spilling breaking and plunging breaking waves propagating over different lee side slopes, it is seen that a steeper slope results in the retention of larger amounts of energy in the fundamental frequency. On the other hand, a milder lee side slope results in lower energy in the fundamental frequency and also lower total energy in the propagating wave. In the context of combining a submerged bar with a floating breakwater, these results indicate that a submerged bar constructed in front of a floating breakwater should have a mild slope sufficient to decompose the incident waves to shorter waves for efficient absorption by the floating breakwater.

Conclusion
The open-source CFD model REEF3D is used to simulate wave propagation over a submerged bar including wave shoaling, breaking and decomposition for regular long waves with T = 2.5 s. The computed free surface elevations at several locations along the length of the flume are compared with experimental data, and generally good agreement is seen both in terms of the phase and the elevation of the free surface variation for both non-breaking and spilling breaking waves. Good representation of wave shoaling and decomposition of the wave on the weather side and lee side slopes, respectively, is obtained in the simulations. The high-order discretization schemes in the model result in realistic modeling of the nonlinear wave interactions and the dispersion characteristics of the decomposing waves for the more challenging case with long waves with T = 2.5 s, showing strong decomposition on the lee side of the bar. Spilling breakers are observed on the bar crest for an incident wave height of 0.042 m (ξ = 1.0623) and on further increase of the incident wave height to 0.052 m (ξ = 0.9547), plunging breakers are observed. A difference of about 24% of the local wave height is seen in the near post-breaking region due to complex air-water interaction. The wave transformation and decomposition is thoroughly analyzed and the following conclusions are made regarding the wave transformation process along the bar.
• Non-breaking waves with higher incident amplitude increase in crest elevation until the end of the bar crest. • Breaking waves with higher incident amplitude increase in crest elevation until the breaking region on the bar crest. • Breaker classification using the surf similarity numbers based on emergent sloping beaches cannot be applied directly in this scenario. • Higher incident waves keep increasing their celerity and gain in wave phase over the lower wave heights until end of the bar crest.
The power spectra of the free surface elevations along the bar provided the following results regarding the wave decomposition process.
• Significant reduction in the wave energy at the fundamental frequency is seen for all the cases simulated and higher incident waves transfer a larger amount of energy to their higher harmonics on the weather side slope. • A distinct pattern is observed in energy transfer amongst the harmonics with the first, second and third harmonics containing their maximum energies at the initial, the final and in the intermediate stages over the bar, respectively. • A steeper lee side slope results in larger energy retention in the fundamental frequency, whereas a milder lee side slope effectively decomposes the incident wave to its second harmonic.
The current study presents an analysis of regular wave propagation over a submerged bar for non-breaking and breaking waves along with the energy and free surface transformations involved. The study can be extended to irregular wave propagation over a submerged bar along with the effect of the weather side and lee side slope of the bar and the length of the bar crest on the wave transformation and decomposition. Further, the interaction of the waves with a combined submerged bar and floating breakwater to provide effective coastal protection with minimum visual impact on the coastline can be investigated.

Disclosure statement
No potential conflict of interest was reported by the authors.