NON-HYDROSTATIC MODELLING OF INFRAGRAVITY WAVES USING SWASH

This paper presents numerical modelling of the near shore transformation of infragravity waves induced by bichromatic wave groups over a horizontal and a slo ping bottom. The non-hydrostatic model SWASH is ass es ed by comparing model predictions with analytical solutio ns over a horizontal bottom and with detailed labor tory observations for a sloping bottom. Good agreement b tween model predictions and data is found througho t the domain for bound infragravity waves. Furthermore th model predicts greater outgoing free infragravity wave-heights for steeper slope regimes which is consistent with the measurements. The model however tends to overes timate the magnitude of the outgoing infragravity waves.


INTRODUCTION
Infragravity (ig) waves are surface gravity waves with periods ranging from 20s to 250s which are generated by the group structure of the short waves.Two main generation mechanisms have been identified: (i) Longuet-Higgins and Stewart (1962) showed the generation of bound ig-waves by shortwave groups through spatial gradients in the radiation stress and (ii) Symonds et al. (1982) demonstrated that oscillations of the breakpoint induced by short-wave groups at the beach generates free ig-waves.
IG-waves are significant for harbour resonance (e.g.Bowers 1977), moored vessel motions (e.g.Naciri et al., 2004), collapse of ice shelves (e.g.Bromirski et al., 2010) and dune erosion (e.g.Roelvink et al., 2009), which makes ig-waves an important subject for coastal engineers.An accurate prediction of ig-wave conditions in a coastal region requires a proper description of the generation and nearshore transformation of ig-waves.In case of complex bathymetries there are no analytical relations to make such predictions and one has to rely on numerical models.
Several types of numerical models are available to study ig-waves.'Surf-beat models' simulate igwaves by combining a wave driver model, which provides the forcing of the wave groups, with a shallow water model that accounts for the near-shore transformation of the ig-waves (e.g.Roelvink et al., 2009).Alternatively, ig-waves can be modelled by phase-resolving models based on a Boussinesqtype formulation (e.g.Madsen et al., 1991), a Reynolds Averaged Navier Stokes (RANS) type formulation (e.g.Lin and Liu, 1998) or a non-hydrostatic approach (e.g.Stelling and Zijlema, 2003).These phase-resolving models account for all relevant near-shore processes (e.g.shoaling, refraction, reflection, and non-linearity) and thereby provide a potentially more accurate, but computationally more expensive approach.
Non-hydrostatic models differ from classical Navier-Stokes models in that the free surface is described by a single valued function that, compared to Navier-Stokes models, allows non-hydrostatic models to efficiently compute free-surface flows.Furthermore their implementation is less complex compared to Boussinesq models thereby improving robustness and maintenance.However, the capability of the non-hydrostatic approach to accurately simulate the generation and transformation of ig-waves has not yet been fully demonstrated, although it has been verified for the nearshore transformation of short-waves (e.g.Ma et al., 2012).We therefore investigated the capabilities of the SWASH non-hydrostatic model (Zijlema et al., 2011) to simulate ig-waves by comparing model predictions with the analytical finite depth equilibrium solution of Longuet-Higgins and Stewart (1960) and with the laboratory observations of Van Noorloos (2003) who investigated the cross-shore transformation of ig-waves over a sloping bottom.
This paper is structured as follows: first the basic equations of the SWASH model are introduced, the relevant details of its numerical implementation are briefly addressed and the boundary condition which was employed in this study is presented.This is followed by the description and results of the two cases considered in this study, (i) bichromatic waves propagating over a horizontal bottom and (ii) bichromatic waves propagating over a sloping bottom.The paper finishes with conclusions.

NUMERICAL MODEL
SWASH (Simulating WAves till SHore) is a hydrodynamic model for simulating non-hydrostatic free-surface flows.It is based on the nonlinear shallow water equations including non-hydrostatic pressure, which are derived from the incompressible Navier-Stokes equations that describe conservation of mass and momentum.In this study we consider unidirectional waves and the equations are therefore presented in a two-dimensional vertical plane.This plane is bounded by the free surface ( , ) , where t is time and x and z are the Cartesian coordinates with z defined upwards and 0 z = located in the still water level.In this framework the equations are given by 0, u w x z where ( , , ) turbulent stresses and ρ is the density of water.An expression for the free surface is obtained by considering the mass balance for the entire water column, which for an incompressible flow results in d 0.
A bottom friction term is added as it is important for ig-waves and during wave run-up.This term prescribes a stress at the bottom boundary which is based on a simple quadratic friction law where f c is a constant friction coefficient, U is the depth averaged current and h d ζ = + is the total water depth.To allow the influence of the bottom friction to extent over the vertical some vertical mixing is introduced by means of a simple turbulent viscosity approximation for the vertical turbulent stresses (e.g.xz v z u τ ν = ∂ , with v ν the vertical eddy viscosity).In this study we employ a constant vertical eddy viscosity of . The numerical implementation is based on an explicit, second-order accurate (in space and time) finite difference method that conserves both mass and momentum at the numerical level.A structured grid is employed to discretize the physical domain.In x-direction the grid has a constant width x ∆ and in the vertical direction the physical domain is divided into a fixed number of layers between the bottom and the free surface.The thickness of the layers is defined in a relative way as a constant part of the water depth, which is similar to -coordinates.A more detailed overview of the governing equations and the numerical implementation is given in e.g.Zijlema et al. (2011).With the numerical implementation used in the SWASH model, good dispersive properties were found even for low vertical resolutions.For only two layers the error in the phase velocity is roughly 1% up to 3 kh ≈ where k is a representative wave number (Zijlema et al., 2011).
The above equations can be directly applied to a quasi-steady breaking bore to estimate overall characteristics such as its energy dissipation.However, this requires a high resolution in both horizontal and vertical direction.To accurately capture wave breaking with a limited amount of vertical layers, we adopt the breaking formulation of Smit et al. (2012).With this formulation the dispersive effects are suppressed in the vicinity of a breaking wave by neglecting the contribution of the nonhydrostatic pressure.A wave is considered to break if t gh , where g is the gravitational acceleration and α is the maximum wave steepness before a wave is considered to break.The absence of dispersive effects ensures that the wave front develops a vertical face and, based on the analogy between a hydraulic jump and a turbulent bore, the model accounts for the energy dissipation by ensuring conservation of mass and momentum.The resulting numerical model, in combination with conditions at the inflow and outlet of the domain, is capable of describing the effects of (nonlinear) shoaling, refraction, diffraction, breaking of waves, wave run-up and nonlinear interactions.

Boundary conditions
Waves are generated at the inflow boundary of the domain where the horizontal velocities are prescribed.In order to include incident bound ig-waves, the velocity signal must be based on higher order wave theory.In case the boundary condition is based on linear wave theory, spurious free waves are generated at the frequency of the bound harmonics.In this study we employ the second-order solution presented by Longuet-Higgins and Stewart (1960), who performed an evaluation of nonlinear wave-wave interactions in finite water depths for two harmonics and proposed an explicit second-order solution for the surface elevation ζ and the velocity potential φ .For bichromatic waves the solution of the surface elevation and velocity potential is given by where subscript 1 and 2 denotes the first-order solution of the individual primary waves, which follows from linear wave theory, and subscript 1,2 denotes the second-order solution due to the non-linear interaction between the two primary waves.The second-order solution is a combination of a sub (i.e.bound ig-wave) and a super harmonic wave, which are forced by the difference and sum interactions, respectively.In this study we focus on the transformation of ig-waves, therefore, we only consider the difference interactions of which the second order solution can be found in the appendix.The incident velocity signal follows from the velocity potential Short-waves are assumed to fully dissipate at a beach whereas infragravity waves can (partly) reflect, therefore, outgoing ig-waves must be absorbed to prevent re-reflections at the boundary.Outgoing ig-waves are detected as the difference between the target surface elevation (6) and the instantaneous surface elevation computed by SWASH.The velocity signal of outgoing ig-waves o u follows from mass conservation in combination with the assumption that outgoing ig-waves are progressive and of constant form where c gd = is the shallow water phase velocity and i ζ is the instantaneous surface elevation.
Outgoing ig-waves are absorbed by adding the velocity due to this motion to the incident velocity signal We employed a radiation condition at the outlet of the domain for simulations with a horizontal bottom.For simulations with a sloping bottom the boundary at the shore is formed by a moving shoreline the (e.g.Zijlema et al., 2011).

BICHROMATIC WAVES OVER A HORIZONTAL BOTTOM
We considered two unidirectional harmonic waves propagating through a domain with a constant bottom depth, where the non-linear interaction between the two primary waves generate a bound wave which is in equilibrium with the short wave forcing.In this manner SWASH is assessed for the generation of bound ig-waves at the boundary and propagation of bound ig-waves through the domain.
Several simulations were performed with constant bottom depths ranging 1-5m.The amplitude of the two primary waves follows from the amplitude depth ratio which is equal for both components and all simulations ( / 0.001 a d = ).The frequency of the primary waves was 1 0.2 Hz 81 f = and 2 0.3 Hz 08 f = for the first and second component, respectively.The range of constant bottom depths results in representative normalized water depths kd ranging 0.6 2 kd < < , where the representative wave number k corresponds to the second primary wave component.
The velocities at the boundary were prescribed with the second-order boundary condition described before, excluding the absorption of outgoing ig-waves.A radiation condition in combination with a sponge layer was employed at the outlet of the domain to minimize the influence of wave reflection.For all simulations the grid size was / 80 , where the wave length λ corresponds to the smallest primary wave length, which resulted in a range of 0.12m 0.33m x < ∆ < .The time step was 0.0125s t ∆ = which corresponds to approximately 250 points per primary wave period.Three vertical layers were employed to obtain good dispersive properties for the kd values encountered.

Method of analysis
The surface elevation signals ( , ) x t ζ were divided into primary-wave and ig-wave components, which are denoted by pw ζ and ig ζ , respectively.The primary-wave signal and ig-wave signal were obtained from a band-pass filter with a frequency resolution 0.0033Hz f ∆ = . The frequency bands were 1 2 1 .9 .1 0 for the primary wave and ig-waves, respectively, where b f is the frequency of the bound wave.These frequency bands contain the energy of the respective components.
To gain insight in the performance of the boundary condition we distinct between bound and spurious free ig-waves.To estimate the energy of the bound and free ig-wave component, the Fast Fourier Transform (FFT) was employed in the spatial domain for 0 . This analysis is based on the difference in wave length between the bound and free component: the length of the bound wave b k is the difference wave number whereas the free wave number f k follows from the difference frequency and the linear dispersion relationship k gd ω = . The energy of an ig-wave component , ig c E was estimated from the wave number spectrum ( ) where hi k and lo k are the higher and lower limit of the wave number band.The wave number limits were chosen based on visual inspection (illustrated in Figure 1 .This method is applied to estimate the bound and free ig-wave energies for the analytical solution and the numerical results.

Results
The model predictions are in good agreement with the analytical solution for both the surfaceelevation and ig-wave surface-elevation (Figure 2).However, there are some discrepancies between the predictions and the analytical solution.A small phase difference between the predicted and analytical primary-wave and ig-wave surface elevation occurs after a distance of approximately three wave groups, as the predicted primary waves propagate slightly faster through the domain.Furthermore, the predicted ig-wave surface elevation shows an oscillation on a larger length scale than the bound igwave, which is associated with the presence of a spurious free wave in SWASH.In Figure 3 a comparison is made for all simulations between (i) the predicted free and predicted bound energy levels and (ii) the predicted and analytical bound energy levels.The results are presented as a function of kd, where the representative wave number k corresponds to the second primary-wave component.Results for an inflow boundary condition based on linear wave theory are included in the graph.The predicted bound energy levels agree with the analytical bound energy levels for a linear and second order inflow condition.Predicted spurious free energies are small for a second order inflow condition compared to the predicted bound energies, in contrast to an inflow condition based on linear wave theory.The second order inflow condition suppresses the generation of spurious free waves.Mechanics Laboratory at Delft University of Technology.The flume is equipped with a piston-type wave board with second-order wave generation and a reflection compensation system to minimize reflections from the wave maker.The bathymetry consisted of an impermeable smooth concrete beach with a constant slope of 1:35.A horizontal bottom part with a still water depth of 0.7m was located between the wave maker x=0m and the toe of the concrete beach x=8.5m (Figure 4).For the experiments the total simulated time was 10min including 5min spin-up time.Surface elevation time series were obtained using 11 wave gauges with a sample rate of 25Hz.Each experiment was run eight times for different wave gauge locations.The measurements were merged into one data set to obtain a high spatial resolution of 0.5m in the shoaling region and 0.3m in the surf zone.

E(k)
Two sets of experiments were conducted, labelled with 'A' and 'B', of which an overview is given in Table 1.In the first set of experiments the primary wave frequencies (f 1 and f 2 ) were varied and this subsequently resulted in a variation of the bound wave frequency (f b ).In the experiments labelled with 'B' the amplitude of the second primary component (a 2 ) and consequently the magnitude of the bound harmonic was varied.
The variation of the bound frequency is associated with varying ig-wave conditions: Battjes et al. ( 2004) found a frequency dependent ig-wave behaviour and related this to the normalized bed slope , where x h is the bed slope, ω is the ig-wave radial frequency and h is a representative depth.Here, similar to Van Dongeren et al. (2007), the representative depth is chosen at the breakpoint.For relative mild slope regimes ( β <0.1) Battjes et al. (2004) found a dominance of incoming bound ig-waves over breakpoint generated ig-waves, observed a large bound-wave amplitude growth in the shoaling zone and found low ig-wave shoreline reflection.In contrast, for relative steep slope regimes ( β >0.45) breakpoint generated ig-waves dominated, incoming bound waves showed a weak amplitude growth and significant ig-wave shoreline reflection occurred.In the present data set the relative bed slope varies between 0.14-0.29 (Table 1).To reproduce the experiments with SWASH, the inflow boundary in the numerical model was located at 6m x = which corresponds to the location of the first wave gauge.The second-order accurate weakly-reflective boundary condition was employed to generate incident bound waves and absorb outgoing free ig-waves.The wave forcing was based on the wave conditions presented in Table 1, however, for all flume experiments a difference between the measured and target amplitude of the first primary wave was observed and therefore this amplitude was set at the measured instead of the target value, which resulted in 1 0.065m a = for all simulations.SWASH is employed with a time step of 0.001s t ∆ = and a grid size of 0.01m, x ∆ = for which a sensitivity analysis for experiment A1 showed sufficient convergence.We employed two vertical layers which is sufficient for the wave dispersion as the representative normalized water depths encountered in the experiments are 1 (10 ) kd O − = .Two remaining parameters, wave breaking parameter α and constant friction coefficient f c , were calibrated for experiment A1.This resulted in 0.01 , which is lower than the value of 0.6 found for random waves by Smit et al. (2012).

Detrended surface elevation signals ( , )
x t were divided into a primary-wave and ig-wave part with a high-pass filter and low-pass filter, respectively, with a cut-off frequency at To gain insight in the incident (bound) and outgoing (free) ig-waves we decomposed the ig-wave signal with the array method, first presented by Battjes et al. (2004) and later improved by Van Dongeren et al. (2007).We used the same number of sensors as Van Dongeren et al. ( 2007) used for the same data set.The rms ig-wave height of the incoming and outgoing ig-wave component follows from ± is the frequency spectrum of the incoming (+) or outgoing (-) ig-wave component, obtained follow from the array method.The frequency limits are 0

Results
First, we present the results of the cross-shore transformation of the rms wave-height for the primary and ig-waves for experiment A1 and A4.These experiments represent the most 'extreme' igwave conditions considered in this study: A1 corresponds to the mildest relative bed slope and A4 corresponds to the steepest relative bed slope.
The cross-shore variation of the primary wave height is predicted accurately for A1 and A4 (Figure 5).Based on these results we distinct between a shoaling region, where minor wave dissipation occurs, and a surf zone, where significant wave dissipation occurs due to the breaking of primary waves.A significant decay in the wave height is observed for x>25m and this is chosen as the transition between the shoaling region and the surf zone.
For the ig-wave height there is a distinct difference between the cross-shore variation of experiment A1 and A4.For A1 the ig-wave height grows towards the shore, with a small oscillation in the magnitude for x<25m.In experiment A4 the ig-wave height has a nodal structure and the overall magnitude increases towards the shore.This nodal structure is associated with occurrence of a standing ig-wave due to the presence of an outgoing ig-wave.For experiment A1 the predicted magnitudes agree with the measurements, although the oscillation in the ig-wave surface elevation is less pronounced.For A4 the predicted nodal structure agrees with the measurements whereas the magnitude of the igwave height is overestimated.Figure 6 shows the cross-shore variation of the rms-wave heights of the incoming and outgoing igwaves for experiment A1 and A4.In both experiments the incoming ig-wave height increases as the bound wave propagates towards the shore whereas the outgoing ig-wave height remains relatively constant throughout the shoaling region as the outgoing ig-wave deshoals towards deeper water.The growth rate of the incident bound wave is greater in experiment A1 than in A4.For A1 the outgoing igwave height is small compared to the incoming ig-wave height, which indicates small shoreline reflection and significant dissipation of ig-wave energy.In contrast, the outgoing ig-wave height is of similar magnitude as the incoming ig-wave height for A4.This is associated with greater shoreline reflection and an increased contribution of breakpoint generated ig-waves for experiment A4.SWASH captures the patterns of the incoming and reflected ig-waves for both experiments.For A1 the magnitudes of the predicted incoming ig-wave heights are in agreement with the measurements throughout the domain whereas the outgoing ig-wave height is overestimated in the surf zone.For A4 the predicted incoming ig-wave height compares well with the measurements whereas the outgoing igwave height is overestimated in the shoaling region.The comparison of the primary-wave and ig-wave heights for all gauge locations and all experiments is summarized in Figure 7, where the predictions versus the measurements are plotted.In the shoaling region most results of the primary waves are located within the 10% error bands whereas 0 0.07 0.14 inside the surf zone greater discrepancies are observed as the primary wave energies are typically over predicted.The results for the ig-waves show greater scatter than the primary waves and a significant amount of the results are located outside the 10% error bands (both in the shoaling region and in the surf zone).These observations are confirmed by the mean relative error (Table 2), which is computed with ( ) Furthermore, a comparison is made between the predicted and measured significant ig-wave heights for the incoming and outgoing component (Figure 8).In the shoaling region and the surf zone most results of the incoming ig-wave heights are located within the 10% error bands.The largest discrepancies occur for greater wave heights for which the incoming ig-wave height is generally underestimated.For the outgoing ig-waves, the results show greater scatter compared to the incoming ig-waves and the outgoing ig-wave height is generally overestimated.To further investigate the discrepancy between the predicted and measured outgoing ig-wave height we make a distinction between the results based on the relative bed slope.The discrepancies between the predictions and measurements in the outgoing ig-wave heights are mainly due to experiment A2,A4 and B1-4.For these experiments SWASH tends to overestimate the outgoing igwave height.The discrepancies observed for experiment A2 and B1-4, which have the same relative bed slope (Table 1), are of similar order of magnitude.Despite of the discrepancies, the results show that the outgoing ig-wave height increases for steeper slope regimes, which is consistent with the frequency dependent behaviour observed in the measurements and reported in literature (e.g.Battjes et al. 2004;Van Dongeren et al.,2007).

CONCLUSIONS
In this study we demonstrated the capabilities of SWASH, a non-hydrostatic phase-resolving wave model, in simulating the transformation of infragravity waves induced by a bichromatic wave-group propagating over a horizontal bottom and over a plane slope.
For bichromatic waves over a horizontal bottom, the model results were compared with the analytical solution presented by Longuet-Higgins and Stewart (1960) for a range of normalized water depths ( 0.6 2 kd < < ).The results showed that SWASH correctly generates bound infragravity waves at the inflow boundary and simulates bound infragravity waves throughout the domain.
Model predictions were compared with a laboratory data set presented by Van Noorloos (2003) for bichromatic waves propagating over a plane slope.Results showed that SWASH is capable of predicting the cross-shore transformation of bound infragravity waves.For the outgoing free infragravity waves, discrepancies between measurements and predictions were observed as SWASH tends to overestimate the magnitude of the outgoing infragravity waves.However, the predicted outgoing infragravity wave heights increased for steeper slopes regimes, which is consistent with the frequency dependent behaviour observed in the measurements and reported in literature.
for the spurious free ig-wave,where k ∆ is the wave number bandwidth which depends on the domain length D

Figure 1 .
Figure 1.The estimation of the bound and spurious free ig-wave energies from the wave number spectrum.The vertical dashed lines indicate the theoretical bound and spurious free wave number, the vertical lines indicate the lower and upper limits of the wave number bands.

Figure 2 .
Figure 2. Spatial variation of the primary-wave surface elevation (upper panel) and ig-wave surface elevation (lower panel) for the simulation with a bottom depth of 3m.Analytical solution (black line); predictions (red markers).

Figure 3 .
Figure 3. Ratio of the predicted free and predicted bound energy Ef,P /Eb,P (circles) and ratio of the predicted and analytical bound energy Eb,P /Eb,A (triangles) for a inflow condition based on linear (light red) and second order (black) wave theory.

Figure 4 .
Figure 4. Bathymetry of the experimental set-up.The horizontal line indicates the still water depth and the vertical lines indicate the locations of the wave gauges.
mean-square (rms) wave height was estimated from the variance of the surface elevation signal indicates time averaging and c ζ refers to either the primary-wave or ig-wave surface elevation.

Figure 5 .
Figure 5. Cross-shore variation of root-mean-square wave height of the primary waves (upper panel) and ig waves (lower panel) for experiment A1 (left panels) and A4 (right panels).Measurements (circles); Model results (red line).

Figure 7 .
Figure 7. Predicted (subscript S) versus measured (subcript M) root-mean-square wave height of the primary waves (left panel) and ig-waves (right panel).Results inside the shoaling region are colored black and results inside the surf zone grey.Solid line, perfect agreement; dashed line, 10% error bands.

Figure 8 .
Figure 8. Predicted (subscript S) versus measured (subcript M) root-mean-square ig-wave height of the incoming (left panel) and outgoing (right panel) components.Results inside the shoaling region are colored black and results inside the surf zone grey.A distinction is made between experiment A1-4 with the following colors: A1 (blue); A2 (red); A3 (green) and A4 (cyan).Solid line, perfect agreement; dashed line 10% error bands.