Performance Assessment of Three Turbulence Models Validated through an Experimental Wave Flume under Different Scenarios of Wave Generation

This study aimed to adjust the turbulence models to the real behavior of the numerical wave flume (NWF) and the future research that will be carried out on it, according to the turbulence model that best adjusts to each particular case study. The k-ε, k-ω and large-eddy simulation (LES) models, using the volume of fluid (VOF) method, were analyzed and compared respectively. The wavemaker theory was followed to faithfully reproduce the waves, which were measured in an experimental wave flume (EWF) and compared with the theory to validate each turbulence model. Besides, reflection was measured with the Mansard and Funke method, which has shown promising results when studying one of the most critical turbulent behaviors in the wave flume, called the breaking of the waves. The free surface displacement obtained with each turbulence model was compared with the recorded signals located at three points of the experimental wave flume, in the time domain of each run, respectively. Finally, the calculated reflection coefficients and the amplitudes of the reflected waves were compared, aiming to have a better understanding of the wave reflection process at the extinction zone. The research showed good agreement between all the experimental signals and the numerical outcomes for all the turbulence models analyzed.


Introduction
Computational modelling is one of the most important tools when working in offshore and coastal engineering applications. Its increasing interest can be seen in academia where an important number of articles have been published in relation to wave generation [1], wave breaking [2], floating wind platforms behavior [3], or wave energy converters performance [4]. More specifically, computational fluid dynamics (CFD) show, supported by the evolution of the hardware, great potential in the study of these topics. This is reflected in the creation of numerical wave flumes (NWFs) [5] and numerical wave tanks (NWTs) [6] in the last years. They are used to support the research in experimental infrastructures but not to substitute them, because CFD applications still cannot replace the experimental work [7].
In most of these studies, the turbulence appears as a vital part of the process. Because of this, a good definition is of a great importance when trying to reduce the total cost of testing prototypes [8]. Thus, the correct selection of the turbulence model is essential in the different methods, such as the Eulerian multiphase volume of fraction (VOF) [9] or the smoothed particle hydrodynamics method methodology followed is defined in Section 3. The turbulence models are presented in Section 4, and then, the results are presented and discussed in Section 5. Finally, in Section 6, some conclusions are presented.

Experimental Wave Flume
All the simulations are replicas of the tests carried out in the EWF of the Energy Engineering Department at the School of Engineering in Bilbao (University of the Basque Country, UPV/EHU). The EWF is 12.5 m long, 0.6 m wide, and 0.7 m high. The walls of the flume are made of laminated and tempered glass and the tank base of stainless steel. The commercial software ASDA-Soft (V5) [32] controls the piston-type wavemaker of the flume by controlling the servo drive and servo motor. The period and amplitude of the paddle are defined in the software, which converts the input data in a sinusoidal curve of progressive acceleration and deceleration of the motor. The motor is connected to a linear actuator, which is attached to the paddle, creating the linear movement of the paddle. Then, the intended wave is generated according to the wavemaker theory [33,34] for a pistontype wavemaker. The wave energy extinction area consists of a self-designed parabolic extinction system, which can be adjusted in height and slope to situate it in the position of maximum wave absorption. The scheme of the extinction system can be seen in Figure 1.  The free surface variation (η (m)) is measured in three different points by using resistive-type probes, the position of which can be modified to adjust them according to the wavelengths of the waves generated. The free surface displacement is recorded by means of an ad hoc LABVIEW in time intervals of 3 ms. It is important to highlight that the probes must be separated in equal distances, as it can be observed in Figure 2, avoiding the half and quarter of the wavelength of study to avoid errors due to their positioning. The first resistive probe is situated at 5.5 m from the paddle to ensure, at least, a minimum distance of two wavelengths of the wave generated. In this study, the wavelengths of all the waves are smaller than 2.75 m, which allow, on the one hand, the waves to be fully developed when reaching the first probe and, on the other hand, to have a maximum flume length for each study. The free surface variation (η (m)) is measured in three different points by using resistive-type probes, the position of which can be modified to adjust them according to the wavelengths of the waves generated. The free surface displacement is recorded by means of an ad hoc LABVIEW in time intervals of 3 ms. It is important to highlight that the probes must be separated in equal distances, as it can be observed in Figure 2, avoiding the half and quarter of the wavelength of study to avoid errors due to their positioning. The first resistive probe is situated at 5.5 m from the paddle to ensure, at least, a minimum distance of two wavelengths of the wave generated. In this study, the wavelengths of all the waves are smaller than 2.75 m, which allow, on the one hand, the waves to be fully developed when reaching the first probe and, on the other hand, to have a maximum flume length for each study.

Numerical Wave Flume
The NWF, unlike the EWF, has a 12-m length, due to the elimination of the volume of water that is behind the wavemaker. As this volume has no influence in the wave generation and reflection, its implementation is unnecessary and would only increase the computational cost of each simulation.

The Mesh
For the two-equation models, a 2-D mesh is generated, while for the LES simulations, a 3-D one is used. In both meshes, the paddle is created as a moving wall in one end of the flume, which follows first-or second-order sinusoidal curves, depending on the wave planned. The condition introduced to the wall in the simulations is a grid velocity, so that the equations must define the velocity of the wall in each time-step. Then, the first-order equation is: where A1 [m] is the amplitude of the principal curve and A2 [m] the amplitude of the secondary one. The meshes can be divided in different areas. The air area (number 1 in Figure 3) has no impact in the wave generation and, therefore, has the biggest cells to reduce as much as possible the computational cost. The area below the extinction system (number 2 in Figure 3) has the same cell size than the air area due to its negligible influence in the experiments. The "deep-water" area is the

Numerical Wave Flume
The NWF, unlike the EWF, has a 12-m length, due to the elimination of the volume of water that is behind the wavemaker. As this volume has no influence in the wave generation and reflection, its implementation is unnecessary and would only increase the computational cost of each simulation.

The Mesh
For the two-equation models, a 2-D mesh is generated, while for the LES simulations, a 3-D one is used. In both meshes, the paddle is created as a moving wall in one end of the flume, which follows first-or second-order sinusoidal curves, depending on the wave planned. The condition introduced to the wall in the simulations is a grid velocity, so that the equations must define the velocity of the wall in each time-step. Then, the first-order equation is: where A 1 [m] is the amplitude of the principal curve and A 2 [m] the amplitude of the secondary one. The meshes can be divided in different areas. The air area (number 1 in Figure 3) has no impact in the wave generation and, therefore, has the biggest cells to reduce as much as possible the computational cost. The area below the extinction system (number 2 in Figure 3) has the same cell size than the air area due to its negligible influence in the experiments. The "deep-water" area is the section below the free surface that always contains water (number 3 in Figure 3). This area has smaller cells than the air area because it is affected by the particles' movement in the waves but is bigger than the free surface and extinction areas also to reduce the computational cost of the simulations. These last areas (number 4 and 5 in Figure 3) have a cell size that depends on the wave generated. The height of the cells is calculated by dividing the wave height (H [mm]) by 20, with this being the biggest cell size that ensures the good definition of the wave. The length of each cell is calculated by having an aspect ratio (AR) of 4, which is the relation between the length and height of the cell: In the walls where an important interaction with the fluid exists, the paddle and extinction system, prism layers are introduced to have a better definition of the interaction. The total height of the free surface area (FS H [mm]) is calculated depending on the wave height too, with the total height of this area being: All the boundaries in the 2-D meshes have a non-slip condition and the bottom and the top have a morphing boundary plane constraint to reproduce a realistic movement of the paddle. In the 3-D mesh, the walls satisfy the non-slip condition except the lateral ones, with a slip condition to delete the friction influence of the walls in the wave generation and propagation. This allows a proper comparison with the experimental flume, in which this friction can be neglected due to its much larger width. section below the free surface that always contains water (number 3 in Figure 3). This area has smaller cells than the air area because it is affected by the particles' movement in the waves but is bigger than the free surface and extinction areas also to reduce the computational cost of the simulations. These last areas (number 4 and 5 in Figure 3) have a cell size that depends on the wave generated. The height of the cells is calculated by dividing the wave height (H [mm]) by 20, with this being the biggest cell size that ensures the good definition of the wave. The length of each cell is calculated by having an aspect ratio (AR) of 4, which is the relation between the length and height of the cell: In the walls where an important interaction with the fluid exists, the paddle and extinction system, prism layers are introduced to have a better definition of the interaction. The total height of the free surface area (FSH [mm]) is calculated depending on the wave height too, with the total height of this area being: All the boundaries in the 2-D meshes have a non-slip condition and the bottom and the top have a morphing boundary plane constraint to reproduce a realistic movement of the paddle. In the 3-D mesh, the walls satisfy the non-slip condition except the lateral ones, with a slip condition to delete the friction influence of the walls in the wave generation and propagation. This allows a proper comparison with the experimental flume, in which this friction can be neglected due to its much larger width.

Reynolds-Averaged Navier-Stokes Equations (RANS)
Every Newtonian and continuum flow can be described using the Navier-Stokes equations. The turbulence is represented by the instantaneous values of the flow properties [8]. The equations for turbulence fluctuation are obtained by Reynolds de-composition, which substitutes instantaneous

Reynolds-Averaged Navier-Stokes Equations (RANS)
Every Newtonian and continuum flow can be described using the Navier-Stokes equations. The turbulence is represented by the instantaneous values of the flow properties [8]. The equations for turbulence fluctuation are obtained by Reynolds de-composition, which substitutes instantaneous velocity and pressure fields for a mean value and a fluctuating component [27]. These components are substituted into the non-dimensional continuity and momentum equations using the time averaging: whereū i is the time-averaged velocity, p is the pressure field, ρ is the density of the effective flow, µ is the viscosity, and -ρu i 'u j ' is the Reynolds operator.

Volume of Fluid Method
The VOF method presented by Hirt and Nichols [9] determines the percentage of the volume of a cell that is occupied by a certain fluid: In this study, two fluids are used, water and air, and the free surface data are determined according to the volume fraction of the water, which is 1 if the cell is full of water and 0 if there is no water in the cell. For those cells where the free surface appears, the volume fraction of the water takes a value between 0 and 1, and this determines the amount of water that is in that cell [9]. The orientation of the free surface is perpendicular to the gradient of the volume fraction that the cell has in relation with its adjoining cells. The physical properties of each fluid are obtained from the fluids and their volume fraction as: where ρ i is the density of each fluid and µ i is their viscosity. The general conservation equation that describes the transport of volume fractions can be defined as: where v is the velocity of the fluid, v g is the grid velocity, and S αi is the source of the ith phase. It is important to highlight that, in our case, the source term is 0.

Aims and Methodology
The aim of the present work was to characterize different scenarios of wave generation to be fully reproduced inside a wave flume, aiming to be able to properly adjust the key parameters of the turbulence models to the real behavior of the numerical models, depending on the turbulence model that best adjusts to each particular case study. The methodology followed in this research can be divided in 4 different stages: first, the experiments were theoretically specified to ensure their possible reproduction by the wavemaker of the EWF. Thereafter, the planned experiments were carried out and the data obtained were processed in MATLAB, and compared with the theoretical values given by the wave theory. Next, the simulations were run to validate the numerical model by comparison with the experimental values. Finally, a comparison of experimental and numerical results when studying reflection accomplished all the objectives of this study.

Experimental Procedures
All the waves were generated following the same methodology. First, the position of the paddle was specified (amplitude) using the ASDA-Soft Software. After this, the period was specified, which determines the time interval needed for each complete movement of the paddle. Since the wave period is the same as the period of the paddle, the theoretical value of the wavelength (λ) was calculated with the dispersion Equation (10), which relates the depth (h), period (T), and wavelength (λ). Then, the velocity propagation of the wave, (c [m/s]), was calculated by (11) in order to determine the time, t E , for which the wave will travel forward and backward along the EWF, until it reaches the position of the last probe (x 3 ) again: Each sensor acquired all the data corresponding to the free surface displacement (η i [m]). When the experiment in the EWF ended, the paddle was situated at 0.5 m from the wall to ensure the same starting point for every experiment.

Data Processing by MATLAB
Once the simulations and the experiments were finished, the data were analyzed by MATLAB to obtain the resulting parameters (H, T, and λ) of the incident wave. The results of the EWF were compared with the theoretical values to ensure the correct generation of the wave. Then, these results were used to calculate the corresponding reflection coefficient.
For the analysis of both the experimental or numerical test, the time domain was divided in three stages, depending on the characteristics of the wave: (i) the initial transitory when the wave is developing, (ii) the time interval in which the wave travels along the flume without any influence of reflected waves, and (iii) when the reflection affects the free surface displacement of the incoming incident wave. Therefore, the time interval corresponding to each stage must be defined. As an example, the t R corresponds to the instant when the reflection phenomenon is first measured on each probe: where L tank is the length of the tank and x i is the position of each sensor. Then, the re-reflection time was calculated with Equation (12) for the position of each probe, respectively. To obtain the parameters of the waves, a specific height condition was created in MATLAB that disregards the data corresponding to the interval in which the wave is still developing. This condition was applied until the height was, at least, 70% of the maximum height before reflection appears. This time value was defined as the start time (t S [s]). In each analysis, the time before t S and after t E were deleted. Then, each run was divided into incident wave and reflection analysis. For the analysis of the incident wave, the data between t S to t R were selected and fitted according to the corresponding wave theory, first (linear, Airy's theory) or second (non-linear, Stokes' theory) order, to obtain the experimental and numerical parameters of that wave (H, T, and λ). The values obtained from this fitting were compared to the theoretical values. It is important to highlight that the fitting avoids the variations from the theory that appear in the free surface displacement data.
Finally, the reflection coefficient (K [-]) was calculated with the method developed by Mansard and Funke [35]. This method uses three wave probes to obtain H at each position, and three phase-shift to calculate the amplitude value of both the incident (A I ) and the reflected wave (A R ). Then, the reflection coefficient is obtained:

Numerical
Based on the experiments theoretically defined, a wave selection was carried out to reproduce numerically waves of different properties for the validation of the models. Then, linear and non-linear waves at the intermediate and deep-water regions were simulated. These simulations were created following the descriptions of Section 2.2. Then, the mesh was converted to 2-D for the k-ε or k-ω models to reduce the computational cost, taking advantage of the 2-D nature of the EWF. For the LES model, STAR CCM+ does not allow a 2-D mesh, so a 5-mm-wide mesh was created to not create an excessive increase in the number of cells. Perpendicular planes were defined to replicate the resistive-type probes of the EWF and measure the free surface displacement at the exact position than in the EWF, for which a user-generated function was imposed to determine the position of the free surface. This function uses the volume fraction of the water that can be defined as: where η sim is the position of the free surface in mm, VOF water is the mean volume fraction of water measure in the study plane, 700 is the total height of the NWF in mm, and h study is the depth of study in mm. Thus, the value of η sim at the beginning of each simulation was 0.

Results Comparison
The results were compared for all the turbulence models and for both the theoretically and experimentally obtained results. This comparison was carried out to study the feasibility of the turbulence models for the wave generation, which was the main objective of this research. In order to minimize the variables, the turbulence models were compared at a constant cell size.

Turbulence Models
To have a better definition of the turbulence in the RANS equations, different models have been created throughout the years. Normally, they are ordered with the number of equations that they used, such as zero-equation, one-equation, or two-equation models. In this article, we used two variations of the two-equation models of k-ε and k-ω, apart from the large-eddy simulations.

The Low Reynolds k-ε Model
The standard k-ε low-Re model has the same coefficients as the standard k-ε model; however, it introduced some damping functions to enable its application in the viscous affected regions near-wall [36], making it more compatible with the law of the wall. However, it does not solve the problem with adverse pressure gradients. In this model, the kinematic eddy viscosity can be defined as: where C µ is a model constant, f µ is a damping function, k is the turbulent kinetic energy, and ε is the turbulence dissipation rate. The mathematical modelling can be written in a boundary layer form: In the standard model, the turbulence kinetic energy is defined as: Turbulence dissipation rate (ε) equation: where the σ k = 1 and σ ε = 1.3 are the Prandtl numbers for k and ε. The damping functions in the standard k-ε low-Re are f 2 : and f µ :

The Shear Stress Transport (SST) Model
First presented by Kolmogorov [37], the most important development of the k-ω model was presented by Wilcox and is the one that the program Star-CCM+ uses as a standard. This model is a two-equation model that has clear advantages in comparison with k-ε, achieving higher accuracies in boundary layers with adverse pressure gradients and in shear and separated flows. This model defines the kinematic viscosity as: where k is the turbulence kinetic energy and ω is the specific dissipation rate. The turbulence kinetic energy (k) is defined as: The specific dissipation rate (ω) equation is: The SST model developed by Menter combines advantages from both k-ε and k-ω. Its complete formulation is in [13], where the author shows that the model works with small variances in the shear-stress and improvements for the predictive accuracy for industrial application, reducing the influence of the user-generated grid.

Large-Eddy Simulation (LES)
In the LES approach, the large-scale turbulence is fully resolved while the sub-grid turbulence is modelled. In comparison with RANS approach, only the small isotropic turbulent scales are modelled. It is extremely useful for applications with high Reynolds numbers. The filtered Navier-Stokes equations of LES can be defined as [38]: However, a subgrid model is necessary to have a good definition. For this, the one developed by Smagorinsky was used. This model uses a mixing-length hypothesis to model the subgrid-scale stresses [36].
It is based on a sufficiently high Reynolds number, which ensures the energy transfer from large to small scales, which are responsible for dissipation. It predicts that eddy viscosity reaches its highest values in areas, for example, solid boundaries, where a strong shear appears. The model provides the mixing-length type formula for the subgrid viscosity: where ρ is the density, ∆ is the length scale, and S is given resolving the velocity field: where C S is the model coefficient and f v is the Van Driest damping function, which in the software can be computed as: where A is another model constant and wall y + is the dimensionless wall distance.

Results and Discussion
Six different waves were used in this study to compare the turbulence models, the defining parameters of which are shown in Table 1. These theoretical values were obtained by combining the Le Méhauté chart [39] and the dispersion Equation (10). The selected wave heights and periods, combined with a wide range of depths in the EWF, from 0.3 to 0.5 m, were selected in order to represent the different waves that appear in the Biscay Marine Energy Platform (BIMEP) [40,41]. The fittings of the free surface displacement signals were done using the incident wave interval, which corresponds to all the data between the t S and t R time values. This fitting, shown in Figure 4, matches almost perfectly with the free surface variation data obtained from the experimental runs. This adjustment constantly achieves a quality of overlapping over 95% in the 3 probes for every wave.  The fittings of the free surface displacement signals were done using the incident wave interval, which corresponds to all the data between the tS and tR time values. This fitting, shown in Figure 4, matches almost perfectly with the free surface variation data obtained from the experimental runs. This adjustment constantly achieves a quality of overlapping over 95% in the 3 probes for every wave. The parameters of the wave were obtained by calculating the average of the parameters of the fittings obtained from the signals of the three probes. Then, the values were compared with the theoretical ones. The results and the corresponding relative errors in relation to the theoretical values are shown in Table 2.  The parameters of the wave were obtained by calculating the average of the parameters of the fittings obtained from the signals of the three probes. Then, the values were compared with the theoretical ones. The results and the corresponding relative errors in relation to the theoretical values are shown in Table 2. Focusing on the relative errors, it can be concluded that the piston-type wavemaker appropriately generates the waves according to the wave theory. Only for the wave 1 were relative errors higher than 6% measured.
The computational turbulence models were then compared, following a similar methodology. In this case, the validation of the models was carried out using the experimental results obtained in the EWF as a reference. This step was essential to compare the quality of the models and its computational cost.
The k-ε turbulence model is the one that was first used in this and previous studies [5]. However, a better definition of the turbulent kinetic energy and the turbulence dissipation rate was obtained to improve the accuracy of the simulations. These values were modified to 1 × 10 −5 for the turbulent kinetic energy and 1 × 10 −4 for the turbulence dissipation rate. These changes help to improve the fitting between experimental and numerical values by defining the crests of the waves better, although troughs are yet to be better defined.
The wave parameters obtained from the simulations using the k-ε turbulence model can be seen in Table 3. A limit of 10% for the wave height and 2% for the period were selected to validate the models. Although the errors are low, always below 2% in period, a better definition of the paddle movement is needed because errors become higher when higher waves are generated. It is important to highlight the satisfactory validation of the equation that the moving wall follows. Next, the wavelength was calculated with the phase-shift of the signals between probes and finally the wave propagation velocity, c. It is worth mentioning that the calculation of the wavelength through the phase-shift introduces a cumulative error from the signal fitting, so it will be further studied for its improvement. Table 3. Results obtained from the simulations carried out using the k-ε turbulence model and relative errors with respect to the experimental values of the incident wave.

H (mm) T (s) λ (m) c (m/s) Error H (%) Error T (%) Error λ (%) Error c (%)
Wave 1  The results obtained for the k-ε model were close to the experimental ones, so that the same mesh was used for the simulations carried out with the SST turbulence model. The value for the turbulent kinetic energy was the same as for the k-ε model, as well as the value of the specific dissipation rate (ω). The similar behavior of both two-equation models reinforces the idea of the use of each model for the problem that best suits each of them. The tendency in the relative errors is the same in both models, as it can be seen in Table 4. However, the SST model has a worst definition in the measurements of the wavelength and celerity, while its definition of the wave height has marginal improvements when comparing it to the k-ε model. As aforementioned, the LES model does not allow simulations in a 2-D mesh to be carried out using Star CCM+. Despite using a 3-D mesh, the same rules of cell size were applied, having the same height and length as the 2-D meshes. The model uses the coefficients C S = 0.1, C t = 3.5, A = 25, and the von Karman constant κ = 0.41. Table 5 shows the values and relative errors of each simulation using this LES model. The adjustments of wave height are worse compared to the two-equation models, but the definition of the wavelength is in the same error interval as the two other models. Because of this, it can be concluded that this model reproduces significantly well the experimentally generated waves. A comparison of the free surface displacement was carried out between experiments and simulations because of the low errors measured. This way allows a direct comparison between incident and reflected waves and observation of whether the extinction system affects the free surface displacement. Figure 5 shows these comparisons between the experimental and the three turbulence models used in the first probe. The fitting is significantly good throughout all the experiment, generating a similar energy dissipation by the extinction system, resulting in a similar free surface displacement at all three points. Figure 5 shows a good definition of the wave crests in all the models, but they fall short on the definition of the troughs, for which the LES models achieve a more constant definition of them.
Once it was demonstrated that the free surface displacement fits properly with the k-ε, SST and LES turbulent models for the incident wave, a statistical comparison of the adjustments of the free surface displacement when reflection appears was developed in order to determine the influence of the energy dissipation process and its turbulence modelling. This comparison focused on the wave height and wave period, as shown in Table 6. From this table, it can be concluded that the k-ε turbulent model obtained promising results, but high error in the height was measured for the wave with the shortest period. Table 6. Comparison of the wave height and period of the resultant wave between the experimental and k-ε model, and the error between them.

H exp (mm) T exp (s) H k-ε (mm) T k-ε (s) Error H (%) Error T (%)
Wave 1  Similar results were obtained with the SST model. While the definition of the wave periods is almost identical, a slight improvement is made in the wave height definition. The results can be seen in Table 7, where wave heights do not surpass a 7% error, while wave periods are always maintained at lower than 2%. This can be due to the improvements introduced by Menter, where the behavior of the adverse pressure gradients is improved and thus, a better definition of the troughs is obtained, reaching, in some cases, a better definition. Once it was demonstrated that the free surface displacement fits properly with the k-ε, SST and LES turbulent models for the incident wave, a statistical comparison of the adjustments of the free surface displacement when reflection appears was developed in order to determine the influence of the energy dissipation process and its turbulence modelling. This comparison focused on the wave  Finally, the LES results are shown in Table 8. Although a better definition of the waves is achieved with this model in deeper depths, it shows limitations when working at smaller depths. However, the tendency of the two other turbulence models when defining the wave period is constant and a better definition of the troughs is reached. Finally, the reflection of the extinction system was studied to obtain more data for a proper comparison of the behavior of the turbulence models. The wave breaking, which is related to reflection, is the most turbulent process seen in the flume and therefore, it can result in higher differences between the turbulent models under investigation. Besides, the reflection coefficients allow comparison of the amplitudes of the incident and reflected waves, which confirms both fittings made to the free surface displacement signal.
The reflection method developed by Mansard and Funke uses both heights and phase-shifts to calculate the amplitudes of the incident and the reflected waves. The reflection coefficient is related to the dissipated energy produced in the extinction system by the turbulent process that occurs. Because of this, the study of the reflection phenomenon has great importance when defining the turbulence model and the NWF. In Figure 6, a comparison of the reflection coefficients measured for all the waves, both numerical and experimental, is shown. Despite the different values measured for the coefficients, it is important to highlight that the errors using this method are magnified by the errors measured in the incident and reflected waves.
A direct comparison of the incident and reflected amplitudes is shown in Figure 7 using the method developed by Mansard and Funke. This provides better results in the study of the amplitude of the incident wave when treating the simulations. This can be related to the high sensitivity of the resistive probes of the EWF, which introduce noise in the signal and makes the measurements in some waves difficult. Further research will be needed to improve the implementation of the method in the experiments. Besides, the small variations in the reflected amplitude are magnified when dividing it by the amplitude of the incident wave. For example, for the wave 1, the incident amplitude of the experimental signal is significantly smaller; however, a marginally smaller amplitude of the reflected wave generates a much smaller reflection coefficient. turbulence model and the NWF. In Figure 6, a comparison of the reflection coefficients measured for all the waves, both numerical and experimental, is shown. Despite the different values measured for the coefficients, it is important to highlight that the errors using this method are magnified by the errors measured in the incident and reflected waves. A direct comparison of the incident and reflected amplitudes is shown in Figure 7 using the method developed by Mansard and Funke. This provides better results in the study of the amplitude of the incident wave when treating the simulations. This can be related to the high sensitivity of the resistive probes of the EWF, which introduce noise in the signal and makes the measurements in some waves difficult. Further research will be needed to improve the implementation of the method in the experiments. Besides, the small variations in the reflected amplitude are magnified when dividing it by the amplitude of the incident wave. For example, for the wave 1, the incident amplitude of the experimental signal is significantly smaller; however, a marginally smaller amplitude of the reflected wave generates a much smaller reflection coefficient.

Conclusions
In this study, three different turbulence models were compared (k-ε, k-ω, and LES) and validated using the EWF of the University of the Basque Country. To the knowledge of the authors, hitherto, this is the first work comparing the three turbulence models to determine which of them reproduces the reflection phenomenon studied in the NWF better. The three models showed good correlations when measuring free surface variation in comparison with the experimental data measured in the laboratory.
A more specific study was made in wave generation, testing if the movement of the wavemaker of the NWF follows the wavemaker theory. The performance of the moving wall showed good results, validating the movement equations that were calculated, and accurately matching the wave period and wave height.
The fitting carried out using MATLAB for the experimental and numerical free surface displacement allowed the parameters of the waves to be obtained and compared with each other, as well as the values theoretically calculated.
The main contribution of this research related to the reflection coefficients that were calculated

Conclusions
In this study, three different turbulence models were compared (k-ε, k-ω, and LES) and validated using the EWF of the University of the Basque Country. To the knowledge of the authors, hitherto, this is the first work comparing the three turbulence models to determine which of them reproduces the reflection phenomenon studied in the NWF better. The three models showed good correlations when measuring free surface variation in comparison with the experimental data measured in the laboratory.
A more specific study was made in wave generation, testing if the movement of the wavemaker of the NWF follows the wavemaker theory. The performance of the moving wall showed good results, validating the movement equations that were calculated, and accurately matching the wave period and wave height.
The fitting carried out using MATLAB for the experimental and numerical free surface displacement allowed the parameters of the waves to be obtained and compared with each other, as well as the values theoretically calculated.
The main contribution of this research related to the reflection coefficients that were calculated to compare properly the three turbulence models between each other, and with the experimentally obtained values. The reflection produced by the extinction system was studied in three different ways. The first one, in which a direct comparison of the free surface variation was made, showed good agreement between the experimental and the numerical signals, despite slight disparities in the troughs. The reflection coefficient study showed promising results, although clear differences were seen in them. The results for the simulations showed strong similarity when comparing the incident wave amplitude measured using the Mansard and Funke method with the value measured with the probes. However, the experimental signals showed significant differences in some waves. This could be due to the noise and will be studied in future articles.
Finally, the three studied turbulence models showed a good performance throughout all the time of the research, although they did not achieve a perfect matching with the experimental results. The aim of the article, the use of three different turbulence models to simulate the behavior of the wave flume, was achieved, creating a strong base for future studies, which will be the focus in the fluid-structure interaction of floating breakwaters or wave energy converters, such as oscillating water columns.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: