The evolution of free and bound waves during dispersive focusing in a numerical and physical flume

Since the introduction of the NewWave theory (Lindgren, 1970), focused wave groups are used in physical and numerical studies to investigate the interaction of marine structures and ships with extreme waves. The propagation of such wave groups is associated with high order nonlinearities that can cause considerable deviations from linear and 2 nd order predictions. Consequently, nonlinear numerical models or laboratory tests are needed to accurately describe the evolution of focused wave groups. In the present study, we validate a widely used two- phase Reynolds Averaged Navier-Stokes (RANS) solver realised in OpenFOAM with experimental results for the propagation of steep focused wave groups, using a newly developed methodology based on the separation of harmonics. This approach allows for accurate focusing of wave groups and in-detail examination of the individual evolution of the high order terms, as well as identifying the source of discrepancies between experiments and numerical models. The wave groups comprise long-crested broadbanded Gaussian spectra of increasing steepness propagating in intermediate water depth. The contribution of the nonlinear harmonics to the crest height and overall shape of the wave are also discussed, together with the effect of nonlinear wave interactions on the free-wave spectrum. The rapid growth of 3 rd and 4 th harmonics near focusing as well as the evolution of the free-wave spectrum, cause departures of up to 29% and 22% from analytic linear and 2 nd order predictions. The present results demonstrate that RANS-VoF solvers constitute accurate models to propagate nearly breaking waves.


Introduction
The accurate definition of a design wave for offshore structures, vessels and coastal structures is vital for their survivability, preventing sea accidents with environmental consequences and human losses (Haver, 2000). For a sea state with a given spectrum, the average shape of the largest and steepest non-breaking wave crests can be represented by a theoretical wave form, which is the normalised autocorrelation function of a random ocean surface based on the underlying spectrum scaled by the crest amplitude (Tucker, 1999). When the distribution of the sea surface elevation follows a Gaussian process, this corresponds to the NewWave model, which has been traditionally used for offshore applications, but building on the deep water results, transient wave groups were also studied for intermediate and shallow water depth . Recently, the validity of NewWave has been confirmed for the coastal zone as well (Whittaker et al., 2016;Whittaker et al., 2017) and the method was used to study wave -seawall interaction problems (Sun and Zhang, 2017). Assuming linearity by omitting nonlinear wave interactions and their effect on the underlying spectrum, a NewWave-type wave form is generated when all the components of a wave group come into phase (Tromans et al., 1991). This has been evidenced by a recent field study where the occurrence of extreme wave crests was found to be linked to the dispersive focusing of the most energetic wave components (Christou and Ewans, 2014). The theoretical background on the NewWave theory and its applications in coastal engineering are explained in more detail in the Appendix.
The majority of the studies regarding the evolution of unidirectional wave groups in experimental and numerical wave tanks (NWTs) demonstrated that dispersive focusing of unidirectional wave groups leads to a wave crest at focus, the shape and elevation of which is not predicted by either linear or 2 nd order wave theory Gibson and Swan, 2007;Johannessen and Swan, 2001;Johannessen and Swan, 2003;Shemer et al., 2007). This is an effect of high order nonlinearities in large transient waves, namely the bound and resonant nonlinearities (Gibson and Swan, 2007). Bound nonlinearities are caused by the emergence of nonlinear harmonics that are phase-locked to the wave group and they tend to sharpen the free surface profile. On the other hand, resonant interactions cause redistribution of energy within the wave spectrum by altering the phases and amplitudes of the linear components of the underlying spectrum and practically new free wave components are generated that satisfy the dispersion relation, as shown by Philips (Phillips, 1960). A noticeable effect of resonant interactions is the deterioration of the quality of focusing for increasing steepness of the wave group, manifested as downshifts of the spatial and temporal focus location Ning et al., 2008). The aforementioned bound and resonant nonlinearities have been shown to result in crest elevations higher than 2 nd order Stokes predictions for unidirectional focused wave groups in deep water (Johannessen and Swan, 2001;Johannessen and Swan, 2003), while for intermediate water depth they yield crest elevations lower than linear predictions (Katsardi and Swan, 2011).
It is worth mentioning that the exact resonant interactions cannot be realised in 1D (unidirectional) propagation, because the resonant conditions of the four-wave interaction, i.e. k 1 þ k 2 ¼ k 3 þ k 4 and ω 1 þ ω 2 ¼ ω 3 þ ω 4 , where k i and ω i are the wavenumber and angular frequency of a wave component respectively, cannot be satisfied (Janssen, 2003). For 1D propagation, non-resonant interactions, which can evolve in short time scales, are of particular importance. Typical examples of such nonlinear effects include the instabilities in regular wave trains reported by Benjamin and Feir (1967) for narrow-banded spectra, also known as BF instabilities. For broadbanded spectra and long-term evolution of 100 T p , where T p is the peak period, resonant interactions tend to increase the bandwidth of the spectrum, as reported by Hasselmann (1962). Complementary to Hasselmann's observations, Gibson and Swan (2007) showed that for focused wave groups, changes to the wave spectrum similar to those of long-term evolution can occur rapidly and locally near the focal location within 3-5 wave cycles.
Consequently, the capacity of a NWT to simulate focused wave groups depends primarily on the accuracy of the numerical dispersion and the accurate calculation of wave-wave interactions. The former directly affects the quality and time of focusing and the latter the shape of the spectrum at focus. A NWT designed for steep waves should be able to account for higher than 2 nd order interactions and if wave breaking is involved, to provide the fully nonlinear solution of the problem. For industrial applications in the oil & gas and offshore renewable energy sectors, 3-Dimensional (3D) RANS solvers combined with surface capturing algorithms, like the Volume of Fluid (VoF) method (Hirt and Nichols, 1981), are the standard for wave-structure interaction problems examined with Computational Fluid Dynamics (CFD) tools. Good alternatives are nonlinear potential flow (NPF) solvers (e.g., Swan, 2003 andNing et al., 2009) and Smoothed Particle Hydrodynamics (SPH) (e.g., Dao et al., 2001), but the latter are not broadly used in industry yet (Lin, 2008). The aforementioned fully nonlinear models provide the solution for the velocity potential and surface elevation, without any division in free and bound waves. On the contrary, weakly nonlinear wave models, solving for example the nonlinear Schr€ odinger equation (NLSE), yield the wave field as the variation of the free-wave regime and the bound waves calculated explicitly from former, but they have inherent limitations for studying fluid-structure interaction problems. Therefore, the use of CFD, SPH and NPF codes has gained ground in industry and research despite the high computational cost. This is especially the case for CFD and SPH, which unlike NPF, can also simulate breaking waves and green water effects.
A widely used and acknowledged open-source realisation of RANS models for industrial application is OpenFOAM, which comprises a CFD package for simulating continuous mechanics problems. Regarding the use of OpenFOAM in coastal and offshore engineering, recent research concerns the modelling of transient wave groups (Bredmose and Jacobsen, 2010;Chen et al., 2014;Higuera et al., 2015), but in most cases discrepancies with experimental results are reported, especially for steep wave groups. In previous studies, the reproduction of the steepest and highest NewWave-type focused wave group was limited mainly by the lack of an appropriate correction methodology for the input signal when a CFD model was used and to a lesser extent by the nonlinearity of the wave model when weakly nonlinear models were employed, e.g., Bateman et al. (2001), Katsardi andSwan (2011), Shemer et al. (2007). Here, we exploit the full capacity of the nonlinear solver to examine the spectral evolution by employing a highly controlled wave generation method  that guarantees accurate focusing of the wave group at a predetermined position in space and time in the NWT and simultaneously by conducting a thorough convergence study. The performance of a RANS/OpenFOAM NWT is compared with experimental measurements for a nearly breaking wave group based on a broadbanded Gaussian amplitude spectrum in intermediate water depth. The measured spectra are decomposed into their linearised part, quadratic sub-and super-harmonics, 3 rd and 4 th order harmonics and their propagation is examined separately. Such a detailed examination of the individual harmonics is useful for practical applications, such as overtopping (Orszaghova et al., 2014) and ringing (Fitzgerald et al., 2014). Wave groups of different steepness are employed to study the relative growth of the nonlinear harmonics.
In the remainder of the paper, the testing conditions and the focusing methodology are presented in Section 2 and the numerical solver is described in Section 3. In Section 4, the validation of the NWT is presented. Also, the spectral evolution and the contributions of high order harmonics to the crest elevation are discussed and compared with analytical solutions. The paper closes with concluding remarks and future work.

Generation of focused wave groups
The most challenging aspect of the accurate generation of focused wave groups is the appropriate selection of the amplitudes and phases of the wave components at the inlet. Previously, linear wave theory (Rapp and Melville, 1990), experimental observations , Zakharov's equation (Shemer et al., 2007) and iterative techniques to calculate the required input phases (Chaplin, 1996) or both the phases and amplitudes (Schmittner et al., 2009) have been proposed. Other advanced methods employed NPF models (Fern andez et al., 2014) and pseudo-third-order corrections (Alford and Maki, 2015). In general, the performance of these approaches reduces considerably as the nonlinearity of the wave group increases.
The new methodology for the highly accurate generation of focused waves  tackles the issues of previous techniques. The main difference with other methods is the use of linearised target spectra instead of the full target spectrum as the initial condition at the wave paddle and the utilisation of spectral/harmonic decomposition. The latter is applied on nonlinear wave records to separate the components of various harmonics corresponding to the orders of Stokes expansion. Johannessen and Swan (2003), among others, combined crest and trough focused waves to extract even and odd harmonics of the signal, but Stagonas et al. (2014) also added positive and negative slope focused waves to clearly separate 1 st (linearised part) from 3 rd and higher order components, and 2 nd sum from 2 nd difference terms, with a similar approach as that used for the harmonics of the forces on cylinders (Fitzgerald et al., 2014). The same approach is adopted here through the following application steps: The desired linearised target spectrum, focus point and time are selected and for the same amplitude spectrum, crest focused (CF), trough focused (TF), and positive and negative slope focused waves are generated.
The linearised part of the spectrum is extracted using a suitable linear combination of the four measured spectra and despite being perturbed by 5 th order terms, it can be easily isolated with filtering. The extracted linearised spectrum is compared with the target spectrum and the amplitudes and phases of the wave components are corrected according to the target values: where α in ; α out ; α tgt are the input, measured and target amplitudes of the components of the spectrum respectively and ϕ in ; ϕ out ; ϕ tgt are input, measured and target phases of the components respectively. The wave components of the corrected linearised spectrum are "propagated backwards" from the phase focus location to the inlet boundary using the linear dispersion relation for the calculation of the phases. The amplitudes of the components are not altered according to linear theory. The procedure is applied iteratively until the amplitudes of the linearised part near the inlet and the phases at focus coincide with the target to the desired accuracy.
The introduction of only the linear components ensures accurate reproduction of the wave group by a linear-wave control signal, even if large nonlinear waves are included in the group. It must be noted that the correction of only the linearised part should be the natural choice, since the 2 nd order bound waves (Sharma and Dean, 1981) and the 3 rd order and higher harmonics are uniquely defined by the free components.
Even though the four-wave decomposition returns the harmonics of the signal, it does not per se reveal their nature, e.g. free, bound and spurious waves. For that, observation of the timeseries of the individual harmonics is required. As demonstrated in Section 4.2, the linear harmonics correspond to free waves and thus, the terms "linear" and "free" waves are used interchangeably hereafter. The comparison of the evolution of the nonlinear harmonics with the evolution of the wave group (Sections 4.3 and 4.4) indicates whether they are bound waves or spurious free waves.
Another considerable challenge for in-detail comparisons between numerical models and experiments, such as the evolution of the freewave spectrum, is to ensure same initial conditions, not affected by the individual characteristics of the wave maker. Here, this is tackled by the use of different locations for the amplitude and phase corrections. Matching the generated amplitude spectrum to the target spectrum relatively close to the inlet and not at the focal location allows the exact target spectrum to be tracked. The measured linearised spectrum at focus is a result of the natural evolution of the target spectrum and thus, it is an appropriate representation of an actual wave group that evolved from a linearised target spectrum. Moreover, if the target spectrum is applied exactly on the inlet, discrepancies may occur, due to inconsistencies among different wave generation methods. Nevertheless, the effectiveness of the methodology can decrease if nonlinear harmonics or reflected waves are included in the corrections, because the former do not propagate according to the linear dispersion relation and the latter propagate in the opposite direction of the examined wave group.
The effectiveness of the correction approach is shown in Fig. 1, where the target and measured amplitude spectrum, phase difference of the wave components and surface elevation at focus are presented for three consecutive iterations (a, b, c). It can been seen that the unwanted spatial and temporal downshifting of the wave group are quickly corrected.
To sum up, the methodology for focusing waves allows for generating identical events in experimental and numerical flumes at a predetermined location. The iterative corrections act as an "auto-correction" of the transfer function of the wavemaker, cancelling its particularities (Buldakov et al., 2017). Therefore, the present approach guarantees that wave groups of the same linearised spectrum can be generated independently in a NWT and experiment and the potential inconsistencies are mostly due to spurious waves, which are different due to different wave generators, but vanish with increasing distance from the wavemaker.

Experimental facility and testing conditions
The experiments were conducted in a 20 m long, 1.3 m deep and 2.5 m wide wave flume at University College London (UCL). The flume has a fixed working depth of 1 m and is equipped with seven flap-type wavemakers, the displacement of which is controlled by a linear transport function to reproduce either a target surface elevation or a target spectrum. A parabolic beach is installed at the other end of the flume to minimise reflections, as depicted in the top panel of Fig. 3. For transient focused wave groups only the reflection of spurious long wave components can contaminate the incoming wave group, but by selecting the focal point at 14.1 m from the wavemaker, the reflected long waves arrive after the wave group has passed. Seven resistive wave gauges (WGs) were used to capture surface elevation in the tank with an accuracy of ±1 mm, the locations of the WGs are shown in Table 2.
The wave conditions were selected to represent a broadbanded Gaussian spectrum, which corresponds to a typical theoretical energy distribution and for which the timeseries of the surface elevation have a compact shape at PF. The selected spectrum has also practical advantages compared to spectra with long high frequency tails, such as JONSWAP. The latter require a cut-off at high-frequency wave components that cannot be effectively generated by the wave paddle and have fast dissipation rates along the flume. For the working depth of 1 m and the selected peak frequency of 0.64 Hz, the generated wave groups travel in intermediate water depth and k p d ¼ 1:75, where k p is the wavenumber of the peak frequency wave component and d is the water depth. This is very close to the value of the well-studied New Year Draupner wave (kd ¼ 1:60) (Walker et al., 2004). In order to examine the effect of increasing nonlinearity in spectral changes and free surface shape at focus, wave groups of different steepness were used, i.e. linear, weakly nonlinear and strongly nonlinear limiting non-breaking, as listed in Table 1, where A Th is the linearly predicted crest amplitude at focus. The steepness is controlled by multiplying the amplitude of the wave components of the spectrum by the same factor, while keeping the same frequency bandwidth and focusing location. The characterization of the nonlinearity used in Table 1 is simply to distinguish the wave groups tested.

Free surface modelling in OpenFOAM
OpenFOAM is an open-source and freely available generalized CFD package for solving continuous mechanics problems. The object-oriented programming used in OpenFOAM allows for modularity of the code and great flexibility to include new libraries, while it facilitates the writing of partial differential equations (PDEs). OpenFOAM provides theoretically unlimited parallelization using the OpenMPI implementation of the message passing interface (MPI) (Open CFD, 2012).
OpenFOAM can solve the 3D Navier-Stokes equations for single or multiphase flows using a RANS approach, comprising the continuity and momentum equations, Equations (2) and (3) respectively (Versteeg et al., 2007). To close the set of equations, a quasi-laminar turbulence model was used, since the waves do not break and no turbulent losses are expected. Consequently, the eddy viscosity, turbulent kinetic energy and the Reynolds Stress are set to zero.
where the bold letters indicate a vector field, ρ is the density of the fluid, U the velocity field in Cartesian coordinates, p Ã is the pseudo-dynamic pressure, g the acceleration of gravity, X the position vector, σ T the surface tension coefficient (0.07 kg=s À2 ), κ c the interface curvature, γ i the fluid phase fraction and μ eff the efficient dynamic viscosity.
To solve the governing equations for incompressible multiphase flows of Newtonian fluids in OpenFOAm, the "interFoam" solver is commonly used. In the present study, interFoam was employed with the PISO (Pressure Implicit with Splitting of Operators) (Issa, 1986) algorithm for the appropriate coupling between the velocity and the pressure.
The free surface is tracked with an advanced two-phase flow technique (Berberovi c et al., 2009) based on the VoF method (Hirt and Nichols, 1981). Each cell of the computational mesh is assigned to a value of the fluid fraction (γ i ), with 0; 1; and 0 < γ i < 1 corresponding to air; liquid and their interface, respectively. VoF has been used for highly distorted free surface flows, like overturning waves (Jacobsen et al., 2012;Bredmose and Jacobsen, 2010), but it can suffer from diffusivity, due to low mesh resolution and the finite thickness of the free surface (Rudman, 1997). OpenFOAM, however, has incorporated an advanced algorithm, MULES (multi-dimensional limiter for explicit solution), that guarantees boundedness of scalar fields, in order to improve the accuracy of the representation of the free surface (Open CFD, 2012).
The time stepping is controlled by the Courant condition (C o ) (Courant et al., 1967), which represents the portion of the cell that the advective flow can cover in one time-step. An additional time controller (alphaC o ) for the interface of multiphase flows is also used to ensure stability. OpenFOAM also includes various numerical schemes for the spacial and temporal discretization of the PDEs, which were selected based on preliminary investigations (Vyzikas et al., 2014), aiming for maximum accuracy and optimal computational cost.

Wave generation and absorption in OpenFOAM
The simulation of waves in CFD requires a special set of boundary conditions that provide the appropriate time-dependent velocity field and surface elevation at the inlet and a non-reflective boundary at the outlet. In OpenFOAM, custom boundary conditions have been built for wave generation and absorption (Morgan et al., 2010;Afshar, 2010), but the most advanced and integrated currently available libraries are waves2Foam (Jacobsen et al., 2012) and IHFOAM (Higuera et al., 2013), both including a wide range of wave theories.
Aiming to reduce the computational cost, all simulations presented here were conducted with a fixed boundary and linear wave generation. Preliminary simulations showed that IHFOAM induces local disturbances at the inlet, which decrease with increasing mesh resolution (Vyzikas et al., 2015) and they are not present in waves2Foam, as seen in Fig. 2a. Nevertheless, an almost identical surface elevation is observed at small distance from the boundary (Fig. 2b).
Regarding wave absorption, radically different approaches are proposed in waves2Foam (Jacobsen et al., 2012) and IHFOAM (Higuera et al., 2013). The first employs relaxation zones within which the solution is partially computed by the governing equations of the computational domain and the target solution on the boundary, ensuring a smooth transition between the fully nonlinear domain and the linearised boundary. The second uses an active wave absorption (AWA) method, for which the fixed outlet boundary acts to cancel the incoming fluxes by taking into account the arriving waves and applying an opposite uniform velocity profile based on Shallow Water Equation (SWE). Due to this assumption, AWA is expected to perform less well for intermediate and deep water depth conditions (Higuera et al., 2013). On the other hand, the efficiency of relaxation zones is proportional to their length, while their use has been associated with an increase of the volume of water in the domain (Jacobsen et al., 2012). Relaxation zones and AWA can also be applied on the inlet boundary to ensure stability. In waves2Foam, it was observed that the generation of very steep waves was impossible in the absence of an even short (5 cm) relaxation zone at the inlet, due to instabilities at the free surface and unphysical water velocities (Vyzikas et al., 2014). Despite the local disturbances at the inlet (Fig. 2a), IHFOAM did not induce similar issues. The final selection between waves2Foam (Jacobsen et al., 2012) and IHFOAM (Higuera et al., 2013) was based mainly on the effective absorption of long waves and the computational cost. Spurious long waves, introduced due to the linear wave generation, precede the wave group and if reflected at the outlet can contaminate the timeseries of the surface elevation. After preliminary tests, it was found that the spurious long waves were better absorbed by the AWA, compared with a relaxation zone of at least 2 Â L p long, where L p is the wave length of the peak frequency wave component (see Table 1). The better performance of the AWA for the long waves can be attributed to the employed SWE approach (Higuera et al., 2013). Regarding the computational efficiency, AWA was found to reduce the computational cost by at least 30%. Accordingly, IHFOAM (Higuera et al., 2013) was preferred for all simulations presented hereafter.

Computational domain
The NWT was designed as a two-dimensional numerical mirror of the physical wave flume at UCL, shown in Fig. 3. The computational mesh consists of three layers: a uniform high resolution middle layer with square cells (aspect ratio, AR ¼ 1) extending 0.2 m below and above the still water level, encapsulating the maximum and minimum surface elevation, and two layers that expand from the middle layer to the bottom (maximum AR ¼ 4) and the top (maximum AR ¼ 2) of the NWT. Square cells are recommended for a highly distorted free surface (Jacobsen et al., 2012), while the grading of the cells results in significant savings of computational resources. The refinement around the interface of the two fluids is common practice in CFD and in OpenFOAM is usually performed with the utility "snappyHexMesh" (Morgan et al., 2010), which does not guarantee always an adequately smooth transition to the coarser part of the domain as the cell grading does (Open CFD, 2012).
The free surface elevation was recorded at 100 Hz with WGs located at identical positions in the numerical and physical domain, listed in Table 2. Twenty two additional WGs were included in the numerical flume for better observing the propagation of the wave group. The main differences between the two flumes were: a) the wave generation, performed by a flap-type moving wave paddle in the laboratory and a stationary boundary in the numerical model and b) the absorption of the waves, achieved by a parabolic beach in the physical flume and an outlet boundary with active absorption in the NWT. These differences can result in inconsistencies between the physical and numerical propagation of certain harmonics of the wave group, as demonstrated in Section 4.

Wave boundary conditions
The NWT is a closed rectangular domain consisting of six walls, each being assigned with appropriate boundary conditions for every variable, namely alpha, which refers to the dimensionless scalar field of the phase fraction γ i in Equation (3), Velocity, which refers to the vector field of the velocity components (m=s) and Pressure, which corresponds to the scalar field of the total pressure minus the hydrostatic pressure (Pa ¼ kg=m=s 2 ). The boundary conditions are listed in Table 3 (Open CFD, 2012). The "empty" boundary conditions of the lateral walls produce no solution for the variables normal to the third dimension, which essentially makes the quasi-3D mesh behave as a 2D mesh. The "IHFOAM" boundary conditions are described below.
The surface elevation at the inlet is calculated according to the corrected spectral shape, which, in the present study, was discretized by 320 equidistant frequency components ranging from 0.0078 Hz to 2.50 Hz. However, since the energy at frequencies higher than 1.5 Hz is zero, a  cut-off frequency was applied in order to reduce the computational cost at the inlet boundary condition using 192 wave components. A large number of wave components for discretising the spectrum guarantees accurate replication of the surface elevation and velocity profile of the wave group at the inlet, which is essential for the accurate propagation of the wave group in the NWT (Ning et al., 2009). The number of components should also be adequately high to minimise overlapping of periodic events by having zero surface elevation before and after the focused wave group. To further reduce the computational effort, only a part of the timeseries between 40 s and 70 s was simulated, for the signal of a repeat period of 128 s and a focusing event at 64 s. For the surface elevation and the velocity profile reconstruction at the inlet, linear superposition was used (Dean and Dalrymple, 1991), as seen in Equation (4) and Equation (5), respectively. Second order wave generation could be used as an alternative in order to minimise the spurious free waves, but linear wave generation was preferred here in order to reduce the computational cost at the inlet boundary and allow for consistent comparisons with the experimental linear wavemaker.
where η is the free surface elevation; u and w the horizontal and vertical velocity components (the normal to the NWT component v ¼ 0); ψ the phase of each wave component i; z the distance from the bottom of the NWT; x ¼ 0 m, horizontal distance from inlet boundary and t the time.

Convergence tests
On achieving grid-independent solution, focused waves were generated for combinations of a R2.5-C0.1, R2.5-C0.2, R5.0-C0.2 and R10-C0.2, where R is the minimum cell size in mm and C is the value of C o , which was selected to be the same as alphaC o . The results of the convergence analysis for the measured surface elevation and extracted harmonics after individual corrections of the linearised harmonics with the methodology  are shown in Fig. 4 for the steepest wave group (see Table 1) at the PF location. Coarser resolutions and higher C o result in significant overestimations of the 3 rd þ higher order harmonic, deeper 2 nd order wave troughs, and an overall increase of the crest amplitude by up to 20%.
The computational time varied from approximately 1, 5, 35 and 50 h for the R10-C0.2, R5.0-C0.2, R2.5-C0.2 and R2.5-C0.1 cases respectively, when the CF strongly nonlinear wave group was examined. Simulations with the same cell size (R2.5) and a larger C o run in about 35% faster, but they are accurate only up to 2 nd order. For the R2.5-C0.1 set-up that was finally employed, the wave length of the peak-frequency component (L p ) was discretized by 1 435 computational cells (L p =dx ¼ 1435). The 20 m long domain consisted of 2.48 M cells and the computational time was 16, 34 and 50 h for 30 s of simulation, depending on the steepness of the wave group. As expected, steeper groups required greater computational effort, due to the higher velocities encountered. All simulations were conducted in parallel on a 16-core Intel Xeon E5-2 650 @ 2.6 GHz, using the simple decomposition method, since the cells were uniformly distributed in the x-direction.

Results and discussion
In this section, the numerical results are compared with the experimental measurements at characteristic locations in the wave flume, first as the total (measured) surface elevation and later as separate harmonics. The harmonics are retrieved with an inverse Fourier transform (IFT) of Table 3 Boundary condition for the NWT in IHFOAM (Open CFD, 2012 the analysed spectrum after the four-wave decomposition. As such, the 1 st , 2 nd order sum and 3 rd order components can be readily isolated, while the separation of 2 nd order difference from 4 th order terms is achieved with high/low pass filtering, since they occupy distinctively different frequency bands (Fitzgerald et al., 2014). The validation of the numerical model hereafter concerns mainly the steepest wave group, as it constitutes the most challenging case with the strongest spectral changes. At the end of the section, summarizing results regarding the contribution of each harmonic to the total surface elevation are shown for the wave groups of different steepness. When relevant, comparisons with linear and second order theory are conducted using both the original (target) and evolved linearised spectrum.

Wave group evolution
During the focusing process, the wave group becomes more compact in shape, characterised by a single steep crest. Then it gradually defocuses, with the longer wave components overtaking the shorter, as seen in Fig. 5d for 29 consecutive locations in the NWT. The WGs are not equidistant, but densely located close to PF. The comparison with the experimental results is performed at three characteristic locations: 1.63 m (AM), 9.4 m and 14.1 m (PF) from the wave paddle, as shown in Fig. 5a, b and c, respectively. The largest differences between the numerical and experimental timeseries, calculated as the absolute difference divided by the local elevation, are observed at 1.63 m, where the numerical wave group appears to have 0.0135 m higher crests and 0.0040 m shallower troughs than the experimental one, corresponding to 13.7% and 5.3% of the maximum and minimum surface elevation at that location, respectively. Numerical predictions remain higher at 9.4 m, but the difference, especially around the main wave crest, reduces drastically to 0.0033 m (2.4%), while the difference at the troughs remains proportionally the same (5.1% or 0.0056 m). At focus (14.1 m), a very good agreement is observed for the main wave crest, with discrepancies of only 0.0032 m or 1.5%. This difference is only 1% of the wave height of the focused wave. Negligible differences are seen for the adjacent troughs (0.3%). It is noteworthy that the high frequency distortions observed at the end of the timeseries at the AM location near the inlet boundary are not present in the time history plots at farther locations downstream.

Evolution of the linear harmonics
The evolution of the 1 st order harmonics follow the linear dispersion relation, if wave-wave interactions are ignored, and thus, it is associated with the propagation of the free waves. As it will be demonstrated in Section 4.5, these components have the highest energy content and also determine the evolution of nonlinear harmonics (Krasitskii, 1994). The agreement between the numerical and physical results is very good at all three locations: AM, 9.4 m and PF, as demonstrated in Fig. 6a, b and c, respectively. Discrepancies reduce gradually, as the wave group focuses, resulting in 0.3% (or 0.0004 m) difference at the main crest and a maximum difference of 4.1% (or 0.0039 m) relative to the local elevation at the troughs at PF. The experimental result is consistently higher than the numerical, apart from the recording near the paddle at 1.63 m.
Compared to the linear wave theory estimation with the target spectrum (see Fig. 6c), the linearised harmonics extracted from the numerical model have a higher crest by about 5% or 0.0081 m, shallower troughs by 13% or 0.014 m and lower secondary crests on either side of the main crest. These differences indicate that changes in the free-wave spectrum have occurred.
As discussed in Introduction, the spectral evolution of the free-wave spectrum in 1D is subject to non-resonant interactions. For broadbanded spectra, as the one examined here, widening of the high frequency side, downshifting of the spectral peak and steepening of the low frequency side are expected (Janssen, 2003;Dysthe et al., 2003). These changes are depicted in Fig. 7a, where measurements from consecutive WGs are given to illustrate the gradual evolution of the spectrum. The vertical line at 0.64 Hz denotes the peak frequency (f p ) of the target spectrum. It is seen that as the wave group focuses, the amplitude of frequencies between 0.64 Hz and 0.92 Hz decreases, while the energy content of frequencies higher than 0.92 Hz increases, extending to frequencies up to 2:5f p . The spectral peak is downshifted to 0.59 Hz, followed by slight energy transfer to frequencies just below f p . These results are in strong qualitative agreement with previous works (Gibson and Swan, 2007;Dysthe et al., 2003;Katsardi and Swan, 2011;Janssen, 2003), confirming that non-resonant interactions change the dispersive properties, i.e. both phase and amplitude, of the free wave components. As a result, the time history of the surface elevation is affected. Previous studies (Gibson and Swan, 2007;Katsardi and Swan, 2011) argued that spectral changes result in higher crest elevation at focus, referred to as amplitude sum of the wave components, only in deep and not in intermediate-shallow water, where the amplitude sum actually decreases. However, the present results do not show a reduction, but a marginal increase of approximately 1% at the elevation of the main crest of the linearised harmonics in intermediate water depth (Fig. 7b). Another effect of the spectral change on the linearised harmonics is a considerable shallowing of the neighbouring troughs by 13% (or 0.0135 m) for the strongly nonlinear group, accompanied with widening and shallowing of the adjacent crests.

Evolution of the 2 nd sum and difference components
The propagation of the 2 nd sum components is presented in Fig. 8d for 16 non-equidistant locations in the NWT, starting at the inlet, together with the scaled wave group, allowing for distinction between the free and bound 2 nd order components. Beyond 5 m from the inlet, part of the 2 nd  T. Vyzikas et al. Coastal Engineering 132 (2018) 95-109 sum components, namely the 2 nd sum free waves, are seen to separate from the wave group and travel slower than the group celerity (c g ). As a result, they do not occur simultaneously with the main group at PF. The agreement between the numerical and experimental results seems to improve gradually from AM to PF, as demonstrated in Fig. 8a-c. At 1.63 m the two timeseries of the 2 nd sum harmonics have important differences of up to 50% and the phases are not in very good agreement. At 9.4 m, the maximum difference drops to 6.8% at the crest and 7.7% at the troughs, with the experimental results being higher. At 14.1 m, the numerical prediction is higher by only 2.6% or 0.0012 m at the crest and has a deeper trough by 6.2% or 0.0021 m. The evolution of a wave group in a nonlinear medium induces 2 nd order bound waves, which can be analytically described by the model of Sharma and Dean (1981) for water waves. Thus, the substantial discrepancies near the inlet (AM) can be attributed to the spurious 2 nd order waves created by the linear wave generation. To confirm that, the analytical 2 nd order solution (Sharma and Dean, 1981) based on the target spectrum is plotted in Fig. 8a for the 2 nd sum harmonics. It is shown that, both in the NWT and in the physical flume, linear wave generation induces spurious effects, which cause considerable deviations  T. Vyzikas et al. Coastal Engineering 132 (2018) 95-109 from theory. Nonetheless, the 2 nd order sum solution based on the target spectrum is in better agreement with the physical and numerical results at PF, thanks to the fact that the spurious waves are not present, because their celerity is smaller than that of the main wave group. Thus, the PF should be selected adequately far from the inlet to allow sufficient propagation of the wave group in the nonlinear domain and separation of the spurious free waves. The 2 nd difference harmonics are associated with a long bound wave, in the form of a set-down under a unidirectional wave, being deepest at PF. At 1.63 m, close to the inlet, Fig. 9a shows that there are significant discrepancies between the numerical model and the experiment, with the latter being in relatively good agreement with the 2 nd order solution based on the target spectrum, but still having half of its amplitude. On the contrary, in the numerical model, a spurious long crest is created, which can only be explained as an artefact of the boundary conditions. In fact, when a moving paddle is used with a first order motion, free parasitic long waves, free displacement long waves and local evanescence modes are generated (Sch€ affer, 1996). In the present study, this is manifested by a spurious free wave crest, which, as seen in Fig. 9d, results in preceding erroneous local elevation. Similar unwanted 2 nd order long waves preceding the wave group in a numerical model was reported by Orszaghova et al. (2014). Towards focusing, at 9.4 m in Fig. 9b, the agreement between the numerical model and the experiment, at least regarding the set-down, improves, with the numerical prediction being deeper by 9.2% or 0.0006 m than the experimental result. The spurious crest is still present, and has the same magnitude as the trough. At PF, the numerical 2 nd difference harmonics are deeper by 13.8% (0.0036 m) than the experimental measurement, with both of them being deeper than the analytic solution based on the target spectrum, as seen Fig. 9c.
A potential issue with these long waves is their ineffective dissipation at the end of the wave flume. As seen in Fig. 9c, a reflected wave appears in the experiment, but only 2 s after the main event, while in the NWT the long wave seems to be effectively absorbed.

Evolution of the 3 rd and 4 th components
The evolution of 3 rd and 4 th order harmonics is similar, as seen in Figs. 10 and 11 and thus it is discussed together. In agreement with the work of Katsardi and Swan (2011), it is confirmed that the contribution of 3 rd order terms is important and additionally, that the 4 th order terms are not negligible, when a steep wave group is examined. Consequently, the two-wave decomposition employed in past studies (Johannessen and Swan, 2003;Gibson and Swan, 2007) would not be adequate here and the four-wave decomposition, which treats the 3 rd order harmonics separately, is more appropriate.
As seen in Figs. 10d and 11d, the 3 rd and 4 th order harmonics consist of free and bound waves, with the former lagging behind the main group after 5 m of propagation and the latter increasing in magnitude only after 12.5 m and more rapidly close to PF. This rapid increase was reported by Johannessen and Swan (2003) and Gibson and Swan (2007). For both higher order terms, strong differences are observed near the inlet (1.63 m) (Figs. 10a and 11a), where the numerical predictions are clearly more energetic, but out of phase with the experimental measurements. After the separation of the spurious waves, smaller discrepancies are observed downstream in the flume. At 9.4 m (Figs. 10b and 11b) a relatively good agreement is achieved for the 3 rd order terms, with maximum differences of 17% at crests and 26% at troughs, but the 4 th order terms are still clearly out of phase. At focus (Figs. 10c and 11c), experimental measurements and numerical predictions are in phase, but the latter are seen to overestimate the elevation of the main crest by 20% (or 0.0039 m) and 33% (or 0.0032 m) for the 3 rd and 4 th harmonic, respectively. Deeper troughs are also observed for the numerical model, but the differences are slightly smaller (16% and 30%).

Summary and comparison with analytical solutions
After presenting the evolution of each harmonic separately for the strongly nonlinear group, in this section, we summarize the results for all the wave groups tested and discuss the contribution of the higher order terms to the final elevation and shape of the wave crest at focus. The results are compared with analytical solutions using the target and the evolved free-wave spectrum shown in Fig. 7. These findings are of practical importance as they summarize the validity of analytical solutions in predicting the shape and the crest elevation of a focused wave group according to its steepness.
First, we discuss the departures from linear and 2 nd order theory for wave groups of different steepness (A Th ¼ 0:050; 0:100; 0:154 m). In Fig. 12a, the measured crest elevation at focus and the extracted Fig. 10. (a, b, c) Comparison of the IFT of the 3 rd order harmonics between the (-) numerical and (⋯) physical results. (d) Comparison of the propagation of CF group (⋯) with the 3 rd order harmonics at various locations (-) and at PF (--) in the NWT. linearised, 2 nd sum and difference harmonics are compared with the theoretical 2 nd order solution, based on the target and the evolved linearised wave spectrum, extracted from the numerical simulation. The first apparent feature is that, in agreement with the experimental observations of , nonlinearity increases with the amplitude of the wave group, as indicated by the disproportional growth of the measured wave crest. Moreover, it is noteworthy that the individual quadratic sub-and super-harmonics extracted from the fully nonlinear solution exhibit small discrepancies compared to the 2 nd order solution   12. (a) Comparison of the measured crest elevation (þ), linearised (∘), 2 nd sum (Â) and 2 nd difference harmonics (□) for the numerical (blue markers), 2 nd order theory with evolved free-wave spectrum (red markers) and 2 nd order theory with original free-wave spectrum (--). (b) Comparison of the relative contribution of each harmonic to the measured crest elevation in the numerical model (blue markers) and experiments (--), linearised (∘), 2 nd sum (Â), 2 nd difference (□), 3 rd (Ã) and 4 th (⋄) order harmonics. (c) Surface profile reconstructed only from linearised and 2 nd order terms in the numerical model and the experiment, compared with 2 nd order theory with evolved and original free-wave spectrum. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 4 Relative difference between the measured surface elevation and analytic solution, calculated as (measured/maximum elevation -1), for wave groups of increasing steepness. based on either the target or evolved spectrum. The results of this graph are presented in Table 4 as (%) difference of the analytical solutions against the fully nonlinear computation. As seen, the theoretical results underestimate the numerical crest elevations by up to 29.3% for the strongly nonlinear group, while the error reduces when 2 nd order theory is used with the evolved free-wave spectrum. The advantage of using the evolved free-wave spectrum becomes apparent for the wave groups of moderate steepness, where for 2 nd order theory the error reduces by six times for the quasi-linear group and two times for the weakly nonlinear group. Consequently, a "2 nd order NewWave" calculation with the evolved free-wave spectrum can provide an accurate estimation of the wave crest elevation for wave groups of moderate steepness and a reasonable description of the simulated wave profile for steep wave groups. Thus, it can be potentially used in combination with other techniques (Walker et al., 2004) in order to estimate higher order contributions. Similar analysis, using the original and locally broadened free-wave spectrum, has been performed in a previous study by Johannessen and Swan (2003), who compared experimental results for unidirectional wave groups approaching breaking in deep water depth with linear and 2 nd order Stokes theory. It was reported that the linear theory prediction of the measured crest elevation was improved from 65% to 81% when the evolved free-wave spectrum was used, while the 2 nd order theory prediction was improved from 75% to 92%. Here, the inclusion of the evolved free-wave spectrum marginally improves the theoretical estimation for the strongly nonlinear group, as seen in Table 4. Such differences between the present results and the study of Johannessen and Swan (2003) may be attributed in physical terms to the differences in wave spectra and water depth used in the two studies, but also to the employed different methods for the separation of harmonics, with the four-wave decomposition guaranteeing here that 3 rd and higher bound waves are excluded from the free-wave spectrum.

Quasilinear
Based on these findings, an attempt to reconstruct the time history of the surface elevation at focus from the extracted harmonics is presented in Fig. 12c. It is observed that, even though 2 nd order theory underestimates the wave crest of the strongly nonlinear group, the use of the evolved linearised wave spectrum improves significantly the prediction of the overall shape of the wave group. In particular, the original free-wave spectrum overestimates the depth of the troughs on either side of the main crest by 0.013 m or 13%, whilst the neighbouring crests are also overestimated and appear to be narrower. On the contrary, the troughs are very well predicted when the evolved free-wave spectrum is used -differences of less than 1% from the numerical calculations are reported -and the predicted adjacent crests nearly coincide with those computed. Fig. 12c also serves as an additional validation of the numerical model up to 2 nd order.
After examining the results to 2 nd order, the 3 rd and 4 th order contributions are included in the analysis. The relative contribution of each harmonic, defined as the maximum measured elevation in the time history of the harmonics over the measured crest elevation at focus (A tot ¼ 0:2178m), is presented in Fig. 12b, confirming that the numerical model overestimates all the high order contributions, as discussed in Section 4.3 and 4.4. The contribution of the linearised wave components to the overall surface elevation decreases with increasing steepness, since more energy is transferred to nonlinear harmonics. The linearised terms account for 95%, 89% and 74% of the total crest elevation for the corresponding wave groups of increasing steepness. Regarding the nonlinear terms, for the quasi-linear group (0.050 m), 2 nd order sum harmonics constitute 7.5% of the measured crest elevation, while 3 rd and 4 th order harmonics are negligible, contributing only 0.8% and 0.2%, respectively. The 2 nd sub-harmonics decrease the crest elevation by 3.7%. On the contrary, when the steepness of the wave group reaches the breaking limit (0.154 m), the 2 nd order sum terms contribute three times more (21.5%) to the crest elevation compared with the quasi-linear group and similarly the 2 nd order difference terms cause a three-times deeper setdown (11.2%). In contrast with, the 3 rd and 4 th order harmonics increase by ten times (8.9%) and twenty times (4.5%), respectively. The latter comparison demonstrates the disproportional increase in nonlinearity of the higher orders.
As a final step, the effect that each harmonic has to the shape of the wave group at focus is examined. As such, in Fig. 13 the wave profile of the strongly nonlinear group is gradually built from its harmonics. Clearly, including high order harmonics leads to a narrower and steeper crest and shallower troughs. In more detail, it is seen that if all the positively contributing terms are considered (linear þ 2 nd sum þ 3 rd þ 4 th order terms) the crest elevation is the highest possible (0.2369 m). If the 4 th order terms are excluded the crest decreases in height to 0.2271 m and its shape is slightly widened. When only the linearised and the 2 nd sum harmonics are considered from the numerical simulation a crest elevation of 0.2078 m is predicted, which is relatively close to the Fig. 13. Gradual reconstruction of the surface profile from linear, linear þ 2 nd sum, linear þ 2 nd sum þ 3 rd order terms, linear þ 2 nd sum þ 3 rd þ 4 th order terms and comparison with measured timeseries (including 2 nd difference terms).
T. Vyzikas et al. Coastal Engineering 132 (2018) 95-109 measured crest (0.2178 m), but considerably wider. An almost identical shape and crest elevation to the actual numerical measurement is predicted when the signal is reconstructed with all up to 4 th order harmonics, indicating that the contribution of 5 th order terms is negligible causing only an additional 0.5% increase of the crest elevation. Moreover, it is shown that the set-down of the 2 nd difference terms almost counteracts the positive effect of the 3 rd þ 4 th order terms and it mainly alters the central crest without influencing the neighbouring troughs noticeably. This counteract justifies the very good overall agreement between experimental and numerical results, despite the overpredictions reported in the latter.

Concluding remarks
This study concerns the propagation of focused wave groups in a NWT designed in the widely used CFD package OpenFOAM, using an active wave generation/absorption method (Higuera et al., 2013). The signal at the wavemaker is corrected using a newly developed methodology , which relies on the decomposition of the fully nonlinear signal to its harmonics and allows for accurate focusing of NewWave-type wave groups. The wave groups form a broadbanded Gaussian spectrum of increasing steepness in intermediate water depth.
First, a thorough comparison with experimental results is performed for a nearly breaking wave group at various locations in the NWT, showing that the agreement improves as the group propagates, and at focus the overall comparison with the experiment is excellent. The evolution of the individual harmonics is examined separately along the flume. The linearised harmonics, which essentially correspond to the free-wave spectrum, are in almost perfect agreement with the experimental results, thanks to the iterative corrections of the methodology. At focus, the nonlinear harmonics are in very good agreement, with discrepancies increasing with the order of the harmonics. The nonlinear harmonics in the numerical model appear to have consistently higher energy, but the 2 nd difference terms counteract the excessive contribution of the 3 rd and 4 th harmonics. Closer to the wavemaker the comparison is less impressive, due to the spurious effects of different wave generation in the experiment and numerical model, which causes different erroneous waves. As the wave group propagates the spurious waves separate from the group and the agreement improves significantly.
Emphasis was put on the evolution of the linearised spectrum, which practically dictates the evolution of the nonlinear harmonics. In accordance with previous studies (Gibson and Swan, 2007;Shemer et al., 2007), the analysis showed gradual broadening of the free-wave spectrum, with energy being transferred mainly to the high frequencies, and a noticeable downshift of the peak frequency. This transformation can be attributed to non-resonant four-wave interactions for unidirectional propagation. In contrast with the gradual evolution of the linearised part, the high order harmonics increase rapidly only when the wave group focuses. The results are also compared with analytical solutions using linear and 2 nd order theory with the original and the evolved free-wave spectrum. The use of 2 nd order theory with the evolved free-wave spectrum can provide a very good estimation of the crest elevation and shape of focused wave groups of moderate steepness, but it underpredicts by approximately 20% the crest height of nearly breaking wave groups. This indicates the considerable contribution of 3 rd and 4 th order harmonics for unidirectional steep focused waves.
Overall, the present work shows that a RANS-VoF NWT in Open-FOAM can propagate very steep waves, provided that the boundary conditions are appropriately modified and the solution is well converged, in contrast with previous studies that reported weaknesses of the VoF method (Higuera et al., 2015). Consequently, from a practical point of view, the achieved accurate propagation of steep wave groups encourages further investigations using RANS-VoF models in combination with the methodology and techniques of the present work, for studying in greater depth the interaction of structures with extreme waves. It was also shown through the convergence analysis conducted that the focusing methodology can provide accurate wave shape at the focus location for computationally more efficient NWTs, which is of high relevance for practical engineering applications. Nevertheless, in order to better establish the use of RANS NWTs for more realistic 3D studies with structures, effort should be made to also increase the computational efficiency of the inlet boundary, for example by exploring the accuracy of alternative boundary conditions (Dimakopoulos et al., 2016) or by discretising the wave spectrum with a smaller number of wave components (Ning et al., 2009). engineering are now discussed.
In spite of the NewWave approach being a valuable tool in wave research, there is no commonly agreed practice for its application in, e.g., offshore and coastal engineering design. One of the main reasons for this is the lack of an established methodology for selecting the amplitude of a wave group modelling an extreme event in a random sea. To this end, let us consider a combined probability of occurrence of a storm with a particular spectrum and duration, and of a wave group with particular amplitude within this storm. As a result, equal probability contours in parametric space can be drawn, connecting parameters of wave groups of equal overall return period. These contours will encompass extreme wave groups within moderate sea states, as well as non-extreme waves of severe states. The higher the level of wave extremity within a given sea state as expressed by the ratio of wave height to significant wave height of the sea state -the better it can be approximated by a wave group provided by the NewWave theory. Waves with lower extremity level, i.e., moderately high waves within severe sea states, would experience higher variability, which can be taken into account by considering variations of the wave groups phase spectra. Wave groups with the same linear amplitude at focus, but different phase spectra, have the same probability of occurrence, but lead to different maximum elevations and potentially, given their shape, to different behaviour and type of interaction with a structure; slope focused waves, for example, may result in a nearly vertical crest slamming against a structure. The accurate generation of such groups is illustrated in Fig. 14, where the amplitude and phase spectra of the linearised components at focus and the fully nonlinear elevations are presented for crest, positive and negative slope focused waves; in particular for phase focusing at 0, ±2π=3, ±π=2, ±π=3, and ±π=6. It is noted, that for the generation of the slope focused waves the corrected input for the crest focused waves is used without any additional iterations of the focusing methodology . At this stage, it should also be noted that for severe sea states the maximum crest elevation (or wave height) is to a certain degree restricted by breaking. In this case, the governing parameter specifying the probability of a wave group is not the actual amplitude of the wave, but the amplitude of its linearised components at focus. The latter may be larger than the former, but in practical applications its upper limit will be empirically set by the inescapable occurrence of breaking. Finally, extreme waves within severe sea states, for example the well-known Draupner wave, can be also well approximated by the NewWave theory (Gibbs and Taylor, 2005). Such waves, however, have very high overall return period that may exceed the return periods used in the current engineering design practices.
Overall, the proposed approach provides control over the event's location and has the capacity to generate short in duration and high in repeatability events in experimental facilities and CFD models. It also allows the production of alternative realisations of events with the same occurrence probability. For practical applications, such as wave-structure interaction problems, the possibility to accurately tune the location of the highest crest with respect to the pre-determined location of the structure is provided for first time. Along the same lines, investigations of the structure's interaction with a variety of wave shapes or of the effects imparted from structural design changes, are also made possible without the need for additional long duration tests, which can be impractical for CFD computations and experimental tests. The target is to develop an extended NewWave theory and apply it to offshore and coastal engineering practice. Emphasising on the accurate reproduction of events selected from long irregular wave time histories, this extended approach is envisaged to further complement currently used experimental and computational techniques.