Introduction

This work is developed as part of the WASANO project, funded by the I-Site NExT initiative, which has its main objective in finding an accurate and efficient description of extreme environmental conditions for ocean engineering. In order to do so, there is the need to define the extreme sea states, as well as be able to reproduce them in ocean water tanks.

In this framework, it appears crucial to have highly reliable and efficient numerical solvers which simulate the propagation of waves in wave tanks. Indeed, they can be integrated with experimental solutions in order to foresee the response of structures at sea. Several approaches exist under the potential flow theory. The Boundary Element Method (see e.g. [START_REF] Grilli | A fully nonlinear model for three-dimensional overturning waves over an arbitrary bottom[END_REF]), which makes use of the Green function to solve boundary integral equations, is widely used. Other approaches involve the discretisation of the whole domain, based on finite difference, like OceanWave3D, or finite elements, like IITM-FNPT (see later). Alternatively, it is possible to solve the problem uniquely in the free-surface, by means of the pseudo-spectral methods, as in HOS-NWT. In this study, an extended comparison of three numerical methods that simulate the propagation of fully nonlinear waves in wave tank under the potential theory is performed. These models have been developped by three of the several WASANO partners and they are referred to as Numerical Wave Tanks (NWT). The first method, developed in Ecole Centrale de Nantes, is based on the Higher Order Spectral model, initially ideated by West et al [START_REF] West | A new numerical method for surface hydrodynamics[END_REF] and Dommermuth and Yue [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF]. The second model, OceanWave3D from Danish Technical University, is dedicated to the simulation of waves using high-order finite difference solver [START_REF] Engsis-Karup | An efficient flexible-order model for 3D nonlinear water waves[END_REF]. Finally, IITM-FNPT, by the Indian Institute of Technology Madras, is considered. The latter is based on a finite element discretisation with triangular elements [START_REF] Siriam | Finite element simulation of nonlinear free surface waves[END_REF]. After having presented the numerical methods behind the three models, the metrics used for the purpose of the study, in terms of implementation and accuracy evaluation, are explained. In particular, the numerical results will be compared to free surface elevation measurements from some experiments performed in the LHEEA wave tank. Two sets of experimental data are available: a first one for less steep waves and a second one for steeper waves with breaking events. The quantification of the error will be done both in the time and frequency domain.

Next, for each case, a convergence analysis is done with respect to the most important parameters, as individuated from the theory. This allows evaluating the range of accuracy reached by the three models. Finally, the efficiency in terms of accuracy and computational effort is evaluated. The optimal configuration is sought for the choice of parameters for each model and, then, compared between each other.

Numerical Methods

General Framework

The considered fluid domain, representing the LHEEA cean wave tank, is rectangular and bidimensional (only unidirectional wave case is studied). The vertical z-axis is oriented upwards, with z = 0 corresponding to the free surface elevation at rest and z = -h (where h is the tank depth) to the bottom. The horizontal direction varies from x = 0 (i.e. the wavemaker position) to x = L x , where L x is the tank length.

In the following, potential flow theory will be assumed. As a consequence, the mass conservation reduces to Laplace equation:

△ϕ = 0 ( 1
)
where ϕ represents the velocity potential, defined as:

⃗ ∇ϕ = ⃗ V (2)
with ⃗ V being the fluid velocity. The system needs to be completed with the boundary conditions, which at the walls of the tank are implemented as free-slip conditions:

∂ϕ ∂⃗ n = 0 (3)
where ⃗ n represents the local normal to the fixed walls, in particular the bottom (z = -d).

On the free surface (z = η), the kinematic and dynamic boundary conditions read:

∂ϕ ∂t = -gη + 1 2 | ⃗ ∇ϕ| ∂η ∂t = ∂ϕ ∂z -⃗ ∇ϕ • ⃗ ∇η (4)
After some developments, the equations can be rewritten in terms of free surface quantities, in particular free surface potential φ(x, z, t) = ϕ(x, η, t). Following [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF], the previous free-surface boundary conditions (Eq (4)) become:

∂ φ ∂t = -gη - 1 2 | ⃗ ∇ φ| + 1 2 (1 + | ⃗ ∇η| 2 )W 2 ∂η ∂t = (1 + | ⃗ ∇η| 2 )W -⃗ ∇ φ • ⃗ ∇η (5)
The quantities of interest for the evolution of the free surface (η and φ) can be, therefore, advanced in time solving the previous equations. In order to do so, one needs to evaluate the vertical velocity at the free surface position W = ∂ϕ ∂z (x, z = η, t). Finally, the boundary condition on the wavemaker wall needs to be accounted for and it is expressed as a no-flux condition at the position of the wavemaker x = X(z, t):

∂X ∂t = ∂ϕ ∂x -⃗ ∇X • ⃗ ∇ϕ (6)

HOS-NWT

As introduced in the previous chapter, computing the vertical velocity at the free surface is a key point to solve the nonlinear FSBC. In HOS-NWT, it is done with the Higher Order Spectral method (HOS) [START_REF] Ducrozet | A modified High-Order Spectral method for wavemaker modeling in a numerical wave tank[END_REF], which consists in a double expansion of the velocity potential ϕ and of the vertical velocity at the free surface W . First, they are both expanded as a power series of η:

ϕ(x, z, t) = ∞ m=1 ϕ (m) (x, z, t) W (x, t) = ∞ m=1 W (m) (x, t) (7)
where the index (m) represents the term with order of magnitude equal to the m -th power of η. Moving forward, a Taylor expansion around the free surface rest position is performed for ϕ (m) , obtaining:

ϕ (m) (x, z = η, t) = ∞ n=1 η n n! ∂ n ϕ (m) ∂z n (x, z = 0, t) (8)
Finally, terms from [START_REF] Ducrozet | A comparative study of two fast nonlinear free-surface water wave models[END_REF] and [START_REF] Ducrozet | A modified High-Order Spectral method for wavemaker modeling in a numerical wave tank[END_REF] are regrouped according to their order of magnitude in η. This results in a triangular system, from which the different orders of the velocity potential and, successively, of the vertical velocity at the free surface are computed. Indeed, one gets:

ϕ (m) (x, 0, t) = - m-1 k=1 η k k! ∂ k ϕ (m-k) ∂z k (x, 0, t) W (m) (x, t) = m-1 k=0 η k k! ∂ k+1 ϕ (m-k) ∂z k+1 (x, 0, t) (9)
HOS-NWT relies on the pseudo-spectral approach. The physical quantities are thus decomposed on eigenfunctions of the computational domain ψ and multiplied by the modal amplitudes. In particular, for a rectangular wave tank one gets:

ψ m (x, z) = cos(k m x) cosh[k m (z + 1)] cosh(k m ) (10)
with k m = mπ/L x . Thanks to this definition of the velocity potential, the Laplace problem (Eq (1) with the bottom boundary condition) is implicitly solved. In order to solve the problem numerically, some discretisations need to be applied. First, the order up to which the vertical velocity is evaluated in ( 9) is finite and called HOS order M . Moreover, the functional Fourier space (defined by the set of wavenumbers that appear in [START_REF] Gao | Numerical simulation of deterministic freak wave sequences and wavestructure interaction[END_REF] is computed up to the finite value N x . Also, in order to avoid the phenomenon of aliasing the zero-padding technique is applied. Therefore, the dealiasing can be performed until a chosen parameter. Another discretisation is required in the vertical direction, since the boundary condition on the wavemaker is applied on N z nodesFinally, the time integration scheme used is the fourth-order Runge-Kutta Cash-Karp with adaptive timestep, for which the time tolerance needs to be defined.

OceanWave3D

Developed by DTU, OceanWave3D is a finite difference method which uses a mixed Eulerian-Lagrangian approach. In a first step, the Laplace equation is solved via the finite difference discretisation, and successively the FSBC (4) are used to advance η and ϕ in time.

Figure 1: σ transformation Following [START_REF] Engsis-Karup | An efficient flexible-order model for 3D nonlinear water waves[END_REF], in order to obtain a time-invariant problem, the following change of coordinate in the vertical direction is applied (see Figure [START_REF] Babanin | Breaking of ocean surface waves[END_REF]):

σ = z + h η(x, t) + h = z + h d(x, t) (11)
where h represents the height of the tank. Therefore, the Laplace problem (defined by ( 1) and ( 2)) becomes:

Φ = φ; σ = 1 ∇ 2 Φ + ∇ 2 σ ∂Φ ∂σ + 2∇σ • ∇ ∂Φ ∂σ + ∇σ • ∇σ + ∂σ ∂z 2 ∂ 2 Φ ∂σ 2 = 0; 0 ≤ σ < 1 ∂σ ∂z + ∇h • ∇σ ∂Φ ∂σ + ∇h • ∇Φ = 0; σ = 0 ( 12
)
where Φ(x, σ, t) = ϕ(x, z, t). It should be noticed that all the spatial derivatives of σ can be directly computed if the free surface elevation is known. The wave elevation η and free surface potential φ being known at the current time-step, the system (12) can be solved in order to obtain Φ. The kinematics of the flow is then evaluated via the chain rule:

⃗ V (x, z) = ⃗ ∇ϕ(x, z) = ⃗ ∇Φ(x, σ) + ⃗ ∇σ ∂Φ ∂z (x, σ) W (x, z) = ∂ϕ ∂z ϕ(x, z) = ∂Φ ∂σ (x, σ) ∂σ ∂z (13)
Similarly to HOS-NWT, when the vertical velocity is computed, one can apply the FSBC (5) and move to the next timestep.

As already mentioned, the Laplace problem ( 12) is numerically solved through a finite difference discretisation. In particular, for our 2D case, the grid is composed of N x point along the free-surface and N σ points in the vertical direction. Choosing r nearby points allows order (r -1) finite difference scheme for the 1D first and second derivatives in (x, σ) to be developed in the standard way using Taylor series expansion at each of the x and σ positions on the grid. More in detail, the stencil is composed of r = α + β + 1 points, where α indicates the number of points in the positive coordinate direction and β in the negative direction from the point of interest. Once the problem is linearised thanks to the finite difference scheme, it is resolved with a left preconditioned GMRES (generalised minimal residual) iterative solution. The vector of the solution is, thus, obtained from an iteration when the residuals reach a preset tolerance. The time integration, a 4-th order Runge-Kutta scheme is used and the influence of the timestep size needs to be evaluated.

IITM-FNPT

Following a similar approach to that explained for OceanWave3D, IITM-FNPT also is based on a Mixed Eulerian-Lagrangian scheme (MEL). The main difference stands in the formulation and resolution of the Laplace problem, here obtained through finite elements method.

In particular, as explained in [START_REF] Siriam | Finite element simulation of nonlinear free surface waves[END_REF], the fluid domain Ω (whose boundaries are denoted by Γ) is divided into finite elements connected with n nodes. Formulating the Laplace equation (1) to the associated boundary conditions leads to the following finite elements system of equations:

Ω ∇N i m j=1 φ j ∇N j dΩ| j,i / ∈Γs = = - Γp N i ∂X ∂t (t) dΓ - Ω ∇N i m j=1 φ j ∇N j dΩ| j∈Γs,i / ∈Γs (14)
where Γ s and Γ p are the free surface and wavemaker boundaries; m is the total number of nodes and the potential inside an element ϕ(x, t, z) can be expressed in terms its nodal potential φ j , as:

ϕ(x, t, z) = n j=1 φ j N j (x, z) (15)
Here, N j is the shape function and n = 3 is the number of nodes in the element. Note that for this code the elements are defined with a triangular shape. Moving on to the next step, the horizontal velocity is calculated by fitting cubic splines to the x coordinate and φ values. Basing on continuity condition, one can write:

δx i 6 ∂ 2 f i-1 ∂x 2 + δx i + δx i+1 3 ∂ 2 f i ∂x 2 + δx i+1 6 ∂ 2 f i+1 ∂x 2 = f i+1 -f i δx i+1 - f i -f i-1 δx i ; i = 2, ..., k -1 (16)
where, in our case, f i represents the potential at the i-th node φ i , f ′′ i its second derivative with respect to x and δx i the horizontal spacing between i-th and (i-1)-th node. The previous equation is solved considering k = 5 and the second derivatives of the potential (in our case f i = φ i ) are obtained. It should also be noted that at the two ends the second derivatives are zero, i.e. the natural spline condition. Always following the cubic spline interpolation, once the second derivatives are known, one can compute the first derivatives as:

2 ∂ 2 f fi ∂x 2 + ∂ 2 f i+1 ∂x 2 = 6 δx i f i+1 -f i δx i - ∂f i ∂x (17)
∂ϕ ∂z | z=η = α 2 -1 φ 1 -α 2 φ 2 + φ 3 α (α -1) (z 1 -z 2 ) (18)
where

α = z 1 -z 3 z 1 -z 2 (19)
Note that if α = 1 the nodes are equidistant and Eq (18) reduces to the standard backwards finite difference. Once also the kinematics of the flow at the current time-step is known, the FSBC conditions (4) are applied to update the quantities of interest. In order to overcome difficulties related to numerical instability and high frequency oscillations, regridding at the free-surface is applied. It consists in redistributing the nodes and the associated values at a variable spacing (fitting a cubic spline) after a selected number of timesteps. The time integration scheme here used is a 4-th order Runge-Kutta. Therefore, other than the grid definition (i.e. horizontal and vertical number of nodes N x and N z ), also the choice of the length of the timestep dt will influence the solution.

Numerical Absorption

An absorbing beach is present in the ocean wave tank in order to decrease the waves reflection. It is, then, necessary to reproduce this phenomenon in the numerical models. For all the three considered model, it is done by adding a pressure damping term in the dynamic FSBC from [START_REF] Bonnefoy | A Fully Spectral 3D Time-Domain Model for Second-Order Simulation of Wavetank Experiments. Part B: Validation; Calibration versus experiments and Sample Application[END_REF]. The damping is selected as zero until x = 0.8L x and then smoothly increased. Following [START_REF] Bonnefoy | A Fully Spectral 3D Time-Domain Model for Second-Order Simulation of Wavetank Experiments. Part B: Validation; Calibration versus experiments and Sample Application[END_REF], this set-up corresponds to the best reproduction of the physical beach of the LHEEA wave tank. Thus, the amount of reflection expected is similar to that obtained in the experiments.

Breaking Models

Due to the assumptions made at the beginning (in particular the potential flow hypothesis or the fact that the wave elevation η is uniquely defined on the vertical coordinate), the considered methods cannot simulate the wave breaking. Therefore, for steeper waves, the need to prevent this phenomenon arises. It can be obtained following different approaches and those implemented in the considered codes will be here presented. It should be noted that at the moment no breaking model exists for IITM-FNPT and, therefore, this method will not be part of the study for breaking case.

HOS-NWT

In HOS-NWT, several breaking models can be applied. The first one here introduced is the so-called Tian model. It is based on the idea of dissipating the energy in order to prevent the breaking from occurring. It is applied when the wave configuration exceeds a predefined criterion. The latter is a kinematic criterion, which compares the water particle velocity at the top of the crest U x to the crest speed C x , following Barthelemy's approach [START_REF] Barthelemy | On a unified breaking onset threshold for gravity waves in deep and intermediate depth water[END_REF]:

B x = U x C x = ∂ϕ ∂x 1 C x > threshold (20)
Note that the threshold is set to 0.85 (see [START_REF] Betsy R Seiffert | Simulation of breaking waves using the high-order spectral method with laboratory experiments: Wavebreaking onset[END_REF]). The criterion is computed at every timestep and if it is verified, a dissipation term is added to the FSBC as an extra viscous term. Following Tian [START_REF] Tian | An eddy viscosity model for twodimensional breaking waves and its validation with laboratory experiments[END_REF], one gets:

∂ φ ∂t = -gη - 1 2 | ⃗ ∇ φ| + 1 2 (1 + | ⃗ ∇η| 2 )W 2 + 2ν eddy ⃗ ∇ • ⃗ ∇ φ ∂η ∂t = (1 + | ⃗ ∇η| 2 )W -⃗ ∇ φ • ⃗ ∇η + 2ν eddy ⃗ ∇ • ⃗ ∇η (21)
with ν eddy is defined as:

ν eddy = β H br L br T br (22)
where H br , L br and T br are computed using:

• k b L br = 24.3S b -1.5 • ω b T br = 18.4S b + 1.4 • k b H br = 0.87R b -0.3 with L b = 2L c , k b = 2π/L b , S b = k b (2H c + H t1 + H t2 )/4 and R b = L 2 /L c (with reference to Figure (3)
, depicting the geometry of the wave just prior to breaking); β is dimensionless and in the following will be set to 0.02. Another model that can be applied is the hyperviscous filter. Its action is based on the assumption of an energy cascade from low to high frequencies, modelling breaking as a turbulent process. Thus, it applies a low-pass filter in the Fourier space. Applying the

∂ φk ∂t = E k -µ k φk ∂η k ∂t = F k -µ k η k (23)
where the subscript k represents the Fourier transform of the considered quantity, E k and F k the Fourier transforms of the right hand side of the FSBC (5) and µ k is given by:

   µ k = ak 2 max |k|-k d kmax-k d 2 if |k| ≥ k d µ k = 0 if |k| < k d . ( 24
)
where k max is given by the spatial discretisation and a and k d are to be set by the user.

For the purpose of this study, k d is chosen as 15k max and a as 2. In case that numerical instabilities occurred for most refined discretisations, a will be increased to 3. Figure [START_REF] Barthelemy | On a unified breaking onset threshold for gravity waves in deep and intermediate depth water[END_REF] shows the evolution of µ of Eq (24) with a equal to 2 and k d = 15k max over an example of wave spectrum. Finally, it also possible to implement a combination of the presented breaking models, which applies the Tian model in the physical space and the hyperviscous filter in the Fourier space.

OceanWave3D

OceanWave3D allows the treatment of breaking events with a local smoothing technique, as explained in [START_REF] Kontos | Numerical Modelling of Deep Water Wave Breaking[END_REF]. Similarly to the Tian model in HOS-NWT, a wave breaking onset is detected with the use of a preset threshold. In this model a dynamic limit is fixed and it is based on the downward acceleration at the free surface:

- ∂W ∂t > γg ( 25
)
where γ is a tuning parameter and g is the gravitational acceleration. According to Babanin [START_REF] Babanin | Breaking of ocean surface waves[END_REF], γ has a theoretical value of 0.5. Where the condition is exceeded, a strong smoothing filter is applied on that point and on the neighboring points to each side. Considering the wave elevation at the i-th horizontal node x i , the new free-surface elevation will be given by:

(η(x i , t)) new = 0.25η(x i-1 , t) + 0.5η(x i , t) + 0.25η(x i+1 , t) (26)
This will result in a change of wave amplitude, as well as a phase shift. Following the theory, in this study γ is fixed to 0.5.

Summary

Experimental Data

The accuracy of the numerical results will be evaluated with respect to a set of experiments performed in the LHEEA ocean wave tank, whose characteristics are shown in Figure [START_REF] Bonnefoy | A Fully Spectral 3D Time-Domain Model for Second-Order Simulation of Wavetank Experiments. Part B: Validation; Calibration versus experiments and Sample Application[END_REF].

The horizontal length is L x equal to 46.4m and the depth is h equal to 5m. It should be noticed that the considered wave tank presents flap type wavemaker. The wave elevation measurements come from two experiments, whose goal is to generate a wavefield with a controlled spectrum. In particular, the aim is to have a JONSWAP type spectrum at a specific target location. Table 3 lists the characteristics of the JONSWAP spectrum, the target location and the probe location. The first experiment (Case 1) does not experience breaking events in the wavefield, while the second one (Case 2) does. Figure ( 6) represents the spectra of the wave elevation signal at the considered probes compared to the target JONSWAP spectrum. 4 Methods Implementation and Evaluation

H s (m) T p (s) γ h(m) λ Target x-coord. (m) Probe x-coord. (m)

Wavemaker Motion Implementation

In order to obtain the numerical results, the boundary condition on the wavemaker (6) needs to be implemented. For all the methods, this condition is only considered at the first order of nonlinearity (i.e. ∂X ∂t = ϕ ∂x ). As already mentioned, the LHEEA wave tank presents flap type wavemakers. Nevertheless, since some of the considered methods can only implement piston type, this configuration will be used for all the simulations for a fair comparison. The experimental wavemaker motions are transformed trying not to vary the wavefield generated. Away from evanescent modes and under linear theory, the wave elevation does not change, but for nonlinear waves some discrepancies arise. From the experiments the harmonics η corresponding to the waves generated are available, given for N finite frequencies. In order to obtain the representation of the wavemaker motion in the frequency domain X, the transfer function of a piston wavemaker (as given, for example, in [START_REF] Gao | Numerical simulation of deterministic freak wave sequences and wavestructure interaction[END_REF]) is applied:

T F = η X = i 4sinh 2 (kh) 2kh + sinh(2kh) (27)
Notice that, similarly to η, also X will be evaluated in for N finite frequencies.

Once the Fourier representation of the wavemaker motion is known, the corresponding signal in the time domain is reconstructed as:

X(t) = N i=0 | X(ω i )| cos ω i t + ̸ X(ω i ) (28)
where ω i corresponds to i-th pulsation.

Errors Evaluation

In this section the metrics used to quantify the difference between the measurements from the experiments and the results of the numerical simulations are presented. In particular, the errors will be computed both in the time domain (with the integrated error and the cross-correlation) and in the frequency domain (with the ISSP error). The integrated error (ϵ 1 ) is based on the absolute difference of the signals in the time domain, which is integrated and evaluated as percentage with respect to the reference wave elevation η ref :

ϵ 1 = |η(t) -η ref (t)|dt |η ref (t)|dt (29)
The difference being computed at every singular point in the time domain, it may lead to significant errors if even a small time shift between the two signals is present.

To avoid this discrepancy, the cross-correlation is used to obtain an estimation of the similarity between the wave elevation signals. Given two signals, the cross-correlation between them is defined as:

C ηη,ref (τ ) = (η ⋆ η ref )(τ ) = ∞ -∞ η(t)η ref (t + τ )dt ( 30
)
where η(t) represents the complex conjugate of η(t) and τ is the lag. The value of the crosscorrelation is proportional to the similarity of the two signals, with the second one being shifted in the time domain by τ . Following [START_REF] Aliyar | Efficiency and accuracy of the domain and functional decomposition strategies for the wave-structure interaction problem[END_REF], the following normalization is then applied:

Cηη,ref (τ ) = 1 C η ref η ref (0)C η ref η ref (0) C ηη,ref (τ ) (31)
Consequently, for a complete similarity the maximum value should be reached for zero lag and it should be equal to 1. Thus, the error is estimated as:

ϵ 2 = 1 -max{ Cηη,ref } (32)
For what is above described, it can be interesting to consider the time lag for which the maximum correlation is found τ max , since it corresponds to a phase shift between the signals.

Finally, for the frequency domain analysis the error is computed from the Fourier transform of the wave elevation signal. For this purpose, the improved surface similarity parameter (ISSP) proposed by Kim et al [START_REF] Kim | Real-time phase-resolved ocean wave prediction in directional wave fields: Enhanced algorithm and experimental validation[END_REF] is computed:

ISSP = |F η (f ) -F η,ref (f )| 2 df [|F η (f ) -Fη,ref (f )| + |F η,ref (f ) -Fη,ref (f )|] 2 df (33)
with F η being the Fourier transform of η and F η,ref that of the experimental wave elevation and F the mean value of F .

5 Accuracy of the Models

Nonbreaking Waves

In order to evaluate the efficiency of each method, a convergence study is performed with respect to the most relevant parameters. In particular, a special focus on the influence of the spatial discretisation is first presented. Afterwards, the influence of the other parameters is also discussed. The accuracy of each simulation is evaluated with the methods presented in section (4.2). Basing on the results that have been obtained and will be presented successively, Figure [START_REF] Ducrozet | A comparative study of two fast nonlinear free-surface water wave models[END_REF] shows a qualitative comparison in the time domain of the three different models with the experimental measurements. Here, the most accurate configuration (in terms of ISSP error) is chosen for each model, in particular:

• HOS-NWT: 50 nodes per peak wavelength, M equal to 6, time tolerance fixed to 10 -9

• OceanWave3D: 50 nodes per peak wavelength, 15 vertical nodes, r fixed to 7, GMRES tolerance 10 -8 , 200 steps per T p

• IITM-FNPT: 50 nodes per peak wavelength, 20 vertical nodes, α fixed to 1.2, 200 steps per T p

HOS-NWT

For HOS-NWT the method is tested for varying HOS order and number of points in the x-direction, corresponding to the number of modes in the Fourier expansion. Regarding the latter, the different values of N x are chosen according to the ratio between the maximum wavenumber k max and the wavenumber corresponding to the predicted peak in the wave spectrum k p . In particular, the analysis is performed for k max equal to (5, 10, 15, 20, 25)k p , which, from the wavemaker input (see Table (3)), one expects being around 1.07rad/m. From the definition of the maximum wavenumber (N x = L x k max /π + 1) this correspond to (10, 20, 30, 40, 50) nodes per peak wavelength. It is to be said that, when trying to move towards higher N x , numerical instabilities occur. The HOS order, instead, is varied from 1 to 7. In this first part, the time tolerance is fixed to 10 -4 and its influence will be discussed later. The effect of the aliasing is not investigated in this study, as full dealiasing will be applied. Moreover, the vertical number of nodes will be set to 9. Figure [START_REF] Ducrozet | A modified High-Order Spectral method for wavemaker modeling in a numerical wave tank[END_REF] shows the results for the ISSP error. The left part of the figure gives the convergence with respect to the horizontal discretisation N x and on the right is the similar study for the HOS order M . As it can be appreciated, the influence of the number of points in the Fourier space is significant, with the error that reaches convergence for k max = 20k p (i.e. 40 nodes per peak wavelength). Regarding the influence of the HOS order, it is found that the error tends to stabilizes after M = 3. This results is linked to the fact that, for lower orders, the four-waves interaction are not considered.

A first conclusion can be drawn by looking at the results for the integrated error ϵ 1 and comparing them to the ISSP error in Figure [START_REF] Ducrozet | A modified High-Order Spectral method for wavemaker modeling in a numerical wave tank[END_REF]. In fact, from Figure [START_REF] Engsis-Karup | An efficient flexible-order model for 3D nonlinear water waves[END_REF] one can see how the the evolution of the error against N x and M is analogous to that of the ISSP error. This is due to the linearity of the Fourier transform, which implies that the numerators in Equation ( 29) and ( 33) are strongly related between each other. In fact, the only difference (which is responsible for the appreciable change of error values) is due to the different normalization. Since this result is confirmed for all the other cases, the error ϵ 1 will be no longer presented. The cross-correlation error (Equation ( 32)) is then computed and results are presented in Figure [START_REF] Gao | Numerical simulation of deterministic freak wave sequences and wavestructure interaction[END_REF] similarly to the case of the ISSP error. It can be noticed that the trend of the two errors is very similar. This is due to the fact that the time lag τ max for which the maximum cross-correlation is found is practically constant and equal to 0s, meaning that there is no appreciable phase shift between the numerical results and the experiments. The only exception to this is the case for M equal to 1 for which τ max results in 0.09s (i.e. 0.045T p ). The variation of the time lag for maximum cross-correlation is shown in Figure [START_REF] Grilli | A fully nonlinear model for three-dimensional overturning waves over an arbitrary bottom[END_REF]. No appreciable change is found with respect to N x .

As already mentioned, the influence of the time tolerance used (i.e. 10 -4 ) should be checked. To do so, N x and M are fixed to a value that is previously found to be converged (40 nodes per peak wavelength and M = 6, respectively) and the ISSP error is then evaluated for different time tolerances. The minimum tolerance considered is 10 -9 , since lower values resulted in being too strict to successfully conclude the simulation.

The different errors computed are represented in Figure [START_REF] Kim | Real-time phase-resolved ocean wave prediction in directional wave fields: Enhanced algorithm and experimental validation[END_REF]. The error seems to decrease for tolerances stricter than 10 -5 . However, it must be precised that the overall variation is negligible (the accuracy improves by around 0.5% between a time tolerance of 10 -2 and 10 -9 ). The optimal choice considering also the computational effort will be investigated later, where the results with tolerance equal to 10 -9 , 10 -7 and 10 -4 are considered.

OceanWave3D

A similar study to that of the previous section is performed for OceanWave3D. Here, three cases are considered (with increasing finite difference order), each one presenting the convergence with respect to the mesh definition. The horizontal number of nodes is progressively increased from 10 to 50 nodes per peak wavelength approximately; while in the vertical direction from 7 to 15 points. Three different cases, depending on the value of r are presented separately: 4-th order, 6-th order and 8-th order of finite difference accuracy. For this first part, the GMRES tolerance is set to 10 -8 and the timestep size is equal to 0.01s (200 steps per T p ), since they both correspond to converged results. Beginning with the 4-th order finite difference (i.e. r equal to 5), the evolution of the error is represented in Figure [START_REF] Kontos | Numerical Modelling of Deep Water Wave Breaking[END_REF] with respect to the horizontal and vertical number of nodes on the left and on the right, respectively. It is noticeable that the error begins to converge for 30 nodes per peak wavelength; while regarding the vertical discretisation the results for 12 nodes can already be considered satisfactory. When looking at the cases for higher finite differences orders (Figure [START_REF] Betsy R Seiffert | Simulation of breaking waves using the high-order spectral method with laboratory experiments: Wavebreaking onset[END_REF] and [START_REF] Siriam | Finite element simulation of nonlinear free surface waves[END_REF] for r equal to 7 and 9, respectively), one can deduce similar conclusions. It should be noticed that for r equal to 9 the results for 7 vertical nodes cannot be computed because the stencil is too large. As expected, increasing the dimension of the stencil r improves the accuracy of the solution (especially from 5 to

), with the effect on the efficiency (considering also the CPU time) that will be discussed in the next section. For a more detailed discussion, Figure [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF] shows the variation of the cross-correlation error over the horizontal and vertical discretisation for 6-th order finite difference scheme. Comparing this trend to that of the ISSP error (Figure ( 14)), one can see that, especially for more refined discretisations, the convergence of the error appear similar.

The only case that presents the greatest phase shift is that with 10 horizontal nodes per λ p , for which τ max results in 0.03s, i.e. 0.01T p (no appreciable variation is found with respect to the vertical discretisation). Besides this, the value of the time lag does not vary greatly, but it progressively decreases until it reaches zero for 50 nodes per peak wavelength, as it can be seen in Figure [START_REF] Tian | An eddy viscosity model for twodimensional breaking waves and its validation with laboratory experiments[END_REF]. Finally, it is also found that, regarding the phase shift, the r parameter stands of minor importance, since its variation does not influence the results. Therefore, the cross-correlation error for the other cases is not presented. As explained, the previous result are obtained using a 4-th order Runge-Kutta scheme with fixed timestep (0.01s, i.e. 200 steps per peak waveperiod) for the time integration. The choice of the timestep size, thus, affects the accuracy of the results, as it is here investigated. Having fixed all the other parameters, Figure (18) shows the evolution of the ISSP error for different number of steps per T p . The convergence seems to occur for timesteps smaller than 0.02s (100 steps per T p ).

Finally, another parameter whose influence needs to be estimated is the GMRES tolerance, which defines the accuracy of the resolution of the Laplace problem. The tolerance is 19), the ISSP error seems to stabilize for tolerances stricter than 10 -4 , ensuring the reliability of the previous results. Nevertheless, it should be noted that the variation in the error is not significant. Similarly to the time tolerance for HOS-NWT, the optimal choice of GMRES tolerance will be discussed in Section (6).

IITM-FNPT

Moving forward, the influence of the mesh definition on the results is also studied for IITM-FNPT. Here, two cases are presented, each one corresponding to a different value of the parameter α which determines the distribution of vertical nodes. In particular, the convergence over the number of horizontal and vertical nodes is first presented for α = 1.2, followed by that for α = 2. In order to better visualize the difference, Figure (20) represents the finite element mesh for the two different cases. As predictable from the definition of the parameter (Eq (19)), a greater value of α corresponds to small vertical spacing for points near the free-surface and large vertical spacing for points away from the free-surface. For both cases, the horizontal discretisation is varied from 10 to 50 nodes per peak wavelength; and the vertical from 9 to 15 nodes. Regarding the time integration, it is performed with a RK4 scheme, with, for this first part, timestep size equal to 0.01s. Moreover, regridding is applied every 10 timesteps, since for higher values numerical instabilities occur. per peak wavelength and 12 vertical nodes. Thus, comparing with the mesh refinement in OceanWave3D, one can deduce that convergence is reached for similar refinement of the grid.

Moving to the case with α = 2, whose ISSP errors are depicted in Figure ( 22), one can appreciate the different behaviour especially with respect to the vertical discretisation, as predictable. Here, in fact, the error does not seem to converge before 20 vertical nodes. Therefore, one can deduce that the choice of α equal to 1.2 results being more accurate in this configuration that represents propagation of waves in deep water. In fact, for this type of waves the nodes distribution should not be too sparse at the bottom, as it is for α equal to 2 (see Figure (20)). The results improve significantly when looking at the crosscorrelation error. The latter is presented in Figure ( 23) only with respect to the horizontal discretisation, as it is found to have the most significant influence, both for the cases of α equal to 1.2 (on the left) and 2 (on the right). The considerable difference in the convergence trend if compared to the ISSP error suggests that a phase shift between the numerical results and the experiments affects the results. However, this influence changes considerably depending on the α parameter. Figure (24) shows the variation of the time lag for the maximum cross-correlation for the different horizontal discretisations, which is found to be the most significant parameter. As it can be seen, after a steep decrease between 10 and 20 nodes per λ p , the phase shift tends to stabilize. Some differences, as previously mentioned, are present between the two values of α. In fact, if α is equal to 1.2, the value of time lag for maximum cross-correlation settles at 0.03s (i.e. 0.02T p ) for the converged result of 50 nodes per peak wavelength and 20 vertical nodes. Instead, for α equal to 2, it is found as 0.09s (i.e. 0.05T p ). Therefore, one may conclude that the improvement on the accuracy given by setting α to 1.2 is also due to a considerably smaller phase shift. 25) represents the evolution of the ISSP error for increasing number of timesteps per peak waveperiod (T p = 1.937s). The length of the timestep is varied from 0.1s to 0.01s, corresponding to 20 and 200 steps per waveperiod, respectively. Similarly to OceanWave3D, It can be appreciated that the error tends to stabilize for the last two cases studied (100 and 200 steps per T p ). Therefore, they will be subject of an efficiency analysis in the section [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF]. It should be noted that, for a fair comparison and to prevent numerical instabilities for larger timesteps, the results here presented are obtained with a regridding process after every step, decreasing the overall accuracy of the method.

Breaking Waves

As for the previous case, in this first part a convergence study is performed with respect to the discretisation in the Fourier domain and the HOS order for HOS-NWT and to the mesh refinement for OceanWave3D. Regarding HOS-NWT, three breaking models are tested: the Tian model only, the hyperviscous filter only and the combination of the two. OceanWave3D, instead, is run using the smoothing technique presented in Section (2.6.2). Similarly to the nonbreaking case, also the influence of the time tolerance for HOS-NWT, GMRES tolerance and timestep for OceanWave3D is investigated. Similarly to the nonbreaking case, a qualitative representation of the waves in the time domain is here displayed. The configurations chosen, since they present the highest accuracy, are:

• HOS-NWT; hyperviscous filter: 50 nodes per peak wavelength, M equal to 5, time tolerance fixed to 10 -9 ; a equal to 2

• HOS-NWT; Tian model: 40 nodes per peak wavelength, M equal to 5, time tolerance fixed to 10 -9

• HOS-NWT; Tian hyperviscous filter: 50 nodes per peak wavelength, M equal to 5, time tolerance fixed to 10 -9 ; a equal to 2

• OceanWave3D; smoothing: 30 nodes per peak wavelength, 15 vertical nodes, r fixed to 5, GMRES tolerance 10 -8 , 200 steps per T p

HOS-NWT

Beginning with hyperviscous filter only, Figure (27) shows the evolution of the ISSP error for increasing number of points in the Fourier space (on the left) and HOS order (on the right). Note that, in this first section, the time tolerance is fixed to 10 -4 . The converging trend can be appreciated for both the parameters. However, one can see that when the value of a is increased to prevent numerical instabilities, the ISSP error also rises. This is due to the stronger dissipation applied by the filter (see Eq (24)). Besides for this, the trends are comparable to those obtained for the nonbreaking case, even if, as expected, the value of the error is overall grater. This similarity is also confirmed by the absence of phase shift between the experimental signal and the numerical solution (equal to 0.01s, i.e. 0.004T p ) for values of M greater than 1, as proved in Figure ( 28), where the evolution of the cross-correlation error can be seen being generally analogous to the ISSP error previously presented. This phase shift is consistent for all the cases tested and for all the breaking models used for HOS-NWT. Thus, for the Tian and Tian hyperviscous models the cross-correlation error will not be presented.

Next, the Tian model is used, with Figure (29) showing the ISSP error for increasing number of points in the Fourier space (on the left) and HOS order (on the right). It should be precised that the Tian model results effective up to 40 nodes per peak wavelength, while for a more refined grid, numerical instabilities arise. While for the HOS order the ISSP error reaches convergence similarly to the nonbreaking case and the hyperviscous filter; regarding the horizontal discretisation, the error does not stabilize. This is due to the above-mentioned upper limit of N x . As expected, the value of the error is greater to that obtained for the nonbreaking case.

Moving forward, the same study is performed with the combination of Tian model and hyperviscous filter. The convergence of the ISSP error is depicted in Figure ( 30) with respect to the number of points in the Fourier space and HOS order on the left and on the right, respectively.

The results are comparable to those obtained with the hyperviscous filter only. It can be noticed, anyway, how here the a parameter is increased to 3 in fewer cases, granting higher accuracy. In particular, this difference occurs for 50 nodes per λ p and fifth HOS order, where the dissipation applied by the Tian model avoids having to use a stronger filter. In any case, it is appreciable how this breaking model generally appears to reach higher accuracy if compared to the hyperviscous filter only. However, a detailed comparison between the two models will be investigated in Section (6). Next, the time tolerance is varied from 10 -9 to 10 -2 in order to study its effect on the error. This is done for the Tian model with hyperviscous filter, having fixed the number of discretisation points to 40 nodes per peak wavelength and the HOS order to 5. Regarding the range of tolerances, for values stricter than 10 -9 numerical instabilities occurred. The results for the ISSP error are shown in Figure (31). The trend looks similar to that found for the nonbreaking case and the same conclusion about its influence can be deduced.

OceanWave3D

The accuracy of OceanWave3D with respect to the mesh refinement is here discussed. The breaking model applied is the smoothing technique described in Section (2.6.2). In this first part, the GMRES tolerance and the timestep are, for the moment, fixed to 10 -8 and 0.01s (220 steps per waveperiod), respectively.

Figure (32) and Figure (33) show the results of the study for 4-th order finite difference scheme (i.e. r equal to 5) for the ISSP and cross-correlation error, respectively. It is to be said that for more accurate finite difference discretisations the smoothing appears to be not sufficient and numerical instabilities occur. A convergence of the ISSP error can be appreciated for the vertical discretisation. Regarding the horizontal discretisation, instead, after a steep decrease of the ISSP error between 10 and 20 nodes per peak wavelength the convergence becomes less evident. In fact, for the highest number of points the error appears to slightly increment. This may be due to the effect of the smoothing technique, since the number of steps at which it is applied increase greatly from 30 nodes per λ p up to higher N x . The more definite converging trend for the cross-correlation error suggests that the application of the breaking model results in a significant phase shift. This is confirmed by the lag for which the maximum correlation is found, shown in Figure (34): 0.01s for 10 nodes per peak wavelength against 0.08s for 50 nodes per peak wavelength (considering 15 vertical nodes). The latter value is doubled if compared to the nonbreaking case, confirming that the phase shift is induced by the smoothing technique. Besides the values, one can appreciate how the trend of τ max is opposite to in the nonbreaking case, progressively increasing for greater values of N x . Moving on, the effect of the timestep size is here discussed. After having fixed the other parameters, the evolution of the ISSP error is computed for different number of steps per waveperiod. Unlike for the nonbreaking case, instabilities arise as the timestep size is greater than 0.03s (around 80 steps per T p ). Figure (35) shows the results.

If compared to the nonbreaking case (Figure ( 18)), it can be seen how a larger number of steps per waveperiod is required to obtain converged results. The results with 110 and 220 steps per T p will be object of an efficiency analysis in Section [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF]. Similarly, the influence of the GMRES tolerance is investigated and the evolution of the ISSP error is depicted in The trend is here similar to that obtained in the nonbreaking case, with the ISSP error that tends to stabilize for tolerances stricter than 10 -3 . Nevertheless, similarly to the nonbreaking case, it should be stated that the total variation of the ISSP error is not significant. The effect on the computational effort is discussed later, where the cases with tolerances equal to 10 -3 and 10 -4 are considered.

Efficiency Comparison

In this section, the efficiency of the different models is investigated in terms of accuracy and computational effort, similarly to what done in [START_REF] Ducrozet | A comparative study of two fast nonlinear free-surface water wave models[END_REF]. For each method, the optimal working configuration is sought. In order to do so, several values of the ISSP errors are fixed and the sets of discretisation parameters which reach them with the smallest CPU time are chosen. Moreover, for OceanWave3D and IITM-FNPT the highest efficiency is looked for also in terms of the parameters r and α, respectively and the timestep size. A similar study is also performed with respect to the GMRES tolerance for OceanWave3D and time tolerance for HOS-NWT.

Nonbreaking Waves

Beginning with the nonbreaking case, the previously presented convergence study is used in order to select the most efficient choices in terms of horizontal discretisation and HOS order for HOS-NWT and mesh refinement for OceanWave3D and IITM-FNPT. In particular, as explained above, different values of ISSP errors are fixed and the set of parameters (in terms of grid refinement) that allows reaching that value with the least computational effort for each case is then sought. Note that, for IITM-FNPT and OceanWave3D, this is done with α equal to 1.2 and r equal to 7, respectively. Tables [START_REF] Barthelemy | On a unified breaking onset threshold for gravity waves in deep and intermediate depth water[END_REF] to [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF] Next, also the best choices for the other parameters previously introduced need to be evaluated. This will be done with a series of comparisons where, for the different configurations of Tables [START_REF] Barthelemy | On a unified breaking onset threshold for gravity waves in deep and intermediate depth water[END_REF] to [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF], the accuracy will be shown with respect to the computational effort. Figure (37) presents the comparison between OceanWave3D with the three different values of r, on the left, and IITM-FNPT with respect to α on the right. Each point in the plot refers to one of the configurations listed in Table [START_REF] Bonnefoy | A Fully Spectral 3D Time-Domain Model for Second-Order Simulation of Wavetank Experiments. Part B: Validation; Calibration versus experiments and Sample Application[END_REF] and [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF]. The CPU time necessary to simulate one peak waveperiod is depicted versus the ISSP accuracy (i.e. the complement to 1 of the ISSP error). It can be seen how the results for 4-th and 9-th order finite differences show similar efficiency up to the maximum accuracy reachable by the first. In any case, 6-th order scheme presents overall higher efficiency and will, therefore, considered for the comparison with the other methods. Regarding IITM-FNPT, it is here evident the improvement especially in terms of computational effort due to setting α to 1.2 (for what previously discussed). Having set the most efficient choices for the stencil dimension and the α parameter, a similar investigation is performed with respect to the timestep size both for OceanWave3D and IITM-FNPT. Again, in terms of mesh refinement, the configuration studied are the same as in the previous analysis, listed in Table [START_REF] Bonnefoy | A Fully Spectral 3D Time-Domain Model for Second-Order Simulation of Wavetank Experiments. Part B: Validation; Calibration versus experiments and Sample Application[END_REF] and [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF]. Figure (38) shows the comparisons for OceanWave3D on the left and IITM-FNPT on the right. For both the methods, it is appreciable how the case with 100 steps per peak period presents a higher efficiency, as the computational effort is significantly decreased but the accuracy is not. Therefore, this configuration will be fixed for the following comparisons. The two last parameters whose effect needs to be studied are the time tolerance for HOS-NWT and GMRES tolerance for OceanWave3D. A similar analysis as for the previous parameters is performed and shown in Regarding HOS-NWT, the slight increase in accuracy does not justify using a time tolerance smaller than 10 -4 , as it implies a significantly higher computational effort. Instead, the two cases studied for the GMRES tolerance do not differ greatly, even if for a tolerance equal to 10 -4 the CPU time is reduced without affecting the accuracy of the results. For both these tolerances, thus, the larger values considered results in being the most efficient (because of their negligible influence on the accuracy as previously discussed) Finally, since the optimal choices with respect to all the significant parameters are chosen, the three models can be compared. This results are shown in Figure (40). It is to be said that here, for HOS-NWT, a maximum timestep size of 0.02s (100 steps per T p ) is set for a more fair comparison with the other models in terms of frequency sampling for the error evaluation. However, it is important to precise that a lower computational effort could be achieved without imposing a maximum timestep size. As it can be seen, IITM-FNPT presents a slightly lower efficiency in terms of accuracy and computational effort. Regarding HOS-NWT and OceanWave3D, the latter is found to require competitive computational effort up to the maximum accuracy it can reach (around 0.75). Nevertheless, HOS-NWT results to reach significantly higher accuracy, approximately 0.77. A more general comparison is depicted in Figure (41), where the CPU time of all the simulations run for this study is plotted in a scatter diagram against the corresponding ISSP accuracy. Similarly as for the optimal configurations, it can be seen how HOS-NWT is able to reach highest accuracy; while if a less precise solution is sought the three methods present closer efficiency.

Breaking Waves

Similarly to the nonbreaking case, the optimal configuration in terms of accuracy and CPU time is sought for each method tested. First, basing on the results of the convergence analysis in Section (5.2), the study is performed for the horizontal discretisation and the HOS order for HOS-NWT and the mesh refinement for OceanWave3D. Different values of the ISSP error are set and the configuration that achieves that with the least computational effort is chosen. This is done separately for the three different breaking models studied (hyperviscous and Tian hyperviscous for HOS-NWT and smoothing for OceanWave3D). The parameters configuration resulting from this analysis is resumed in Tables [START_REF] Ducrozet | A comparative study of two fast nonlinear free-surface water wave models[END_REF] to [START_REF] Gao | Numerical simulation of deterministic freak wave sequences and wavestructure interaction[END_REF] Similarly, for the GMRES tolerance in OceanWave3D, increasing the value to 10 -3 improves the required CPU time, without affecting the final error and it will, then, be considered in the following study.

Having set the optimal choice for the parameters above treated, it is possible to visualize a comparison of all the different methods and breaking models, as shown in Figure (44). Like for the nonbreaking case, it should be noted that a maximum timestep size of 0.01s (220 steps per T p ) is set for HOS-NWT for a fair comparison with OceanWave3D, but a smaller CPU time could be achieved without this limitation. It is appreciable how OceanWave3D with the smoothing technique presents good results in terms of computational effort, but considerably smaller accuracy than the other models. This may be related to the limited dimension of the stencil for the finite difference discretisation achievable or to the phase shift induced by the breaking model, as discussed in Section (5.2.2). Regarding the different breaking models tested for HOS-NWT, one can see that they all present similar computational effort up to 0.75 ISSP accuracy. However, the combination of the Tian model and the hyperviscous filter results in reaching a greater value of maximum accuracy. the CPU time of all the cases tested (not only the optimal configurations) is plotted versus the ISSP accuracy.

Conclusion

In this study, an extensive comparison between the three models is presented. First, the theoretical framework, including the numerical methods used for each models is detailed, with a special focus on the techniques to prevent breaking events. After having discussed the metrics used to compare the numerical results with the experiments, a convergence study is performed with respect to the most relevant parameters for each model. Regarding the spatial discretisation (number of points in Fourier space and HOS order for HOS-NWT and mesh refinement for OceanWave3D and IITM-FNPT), the convergence of the errors is used to draw several conclusions. First, it is found that the ISSP error (Eq (33)) presents the same information of the integrated error expressed as in Equation ( 29), but with an improvement related to the normalization. Moving on, regarding the nonbreaking case, for HOS-NWT the number of horizontal discretisation points stand of great importance; while for the HOS order the error stabilizes already after M equal to 3. This is related to the fact that for smaller orders, the four waves interactions are not taken into account, affecting significantly the accuracy of the result. OceanWave3D presents an evident convergence over both the horizontal and vertical number of nodes for all the cases tested. The influence of the order of accuracy of the finite difference scheme is also studied by looking at the evolution of the error with respect to the mesh refinement, and is found to present a considerable improvement in the solution accuracy between the 4-th and 6-th order. For IITM-FNPT, the difference between the two cases studied for different values of α is considerable. For α equal to 1.2, the convergence of the ISSP error is more clear and the error stabilizes for the same mesh refinement as OceanWave3D. The results for α equal to 2 (less vertical spacing for nodes close to the free-surface) show that the ISSP error has issues in reaching convergence. The study of the cross-correlation error presents interesting results regarding the presence of a phase shift between the numerical solution and the experiments. In fact, if the latter is not found for HOS-NWT (except for M = 1), a small value of time lag for the maximum correlation results for OceanWave3D, even if it does not affect greatly the accuracy of the solution, since it progressively decreases for more refined grids. Instead, in IITM-FNPT, with α set to 2, τ max has a considerably greater values than in the other models, indicating a phase shift that affects the accuracy of the solution. For the other parameters, the timestep size results in having a considerable influence on the ISSP error both for OceanWave3D and IITM-FNPT; while the time tolerance in HOS-NWT and GMRES tolerance in OceanWave3D appear to be less relevant. Regarding the breaking case, relevant conclusions can also be deduced. First, for the models available in HOS-NWT, the Tian model is slightly less efficient in terms of stability. The hyperviscous filter is able to reach higher values of N x and M , even if in some cases the filtering needs to be strengthened by increasing a and, thus, increasing the error. The combination of the two above-mentioned models appears to have a better stability and a clearer convergence trend. No appreciable phase shift results for any of these models. The same cannot be said for the smoothing technique of OceanWave3D, where a significant difference exists between the evolution of the ISSP and cross-correlation error over the mesh refinement. Indeed, the time lag for the maximum correlation presents an opposite trend to that of the nonbreaking case, increasing for high N x , when the smoothing is more consistently applied. Also, the dissipation induced by the smoothing is not sufficient for finite difference accuracy orders greater than 4. It is also found that smaller timestep sizes are required if compared to the nonbreaking case.

Instead, not a great difference results for the convergence over the time tolerance (HOS-NWT) or the GMRES tolerance (OceanWave3D). It should be precised that, in this work, the influence of some parameters has not been taken in consideration. For instance, the number of steps after which the regridding is applied in IITM-FNPT is here set as great as possible until numerical instabilities occur. Similarly, the dealiasing parameter is always taken as for full dealiasing.

Next, the computational effort is taken also into account and the optimal configuration for each model is sought. The study is performed considering the selected configurations of spatial discretisation basing on their efficiency for each model. Each parameter is, then, studied separately. In this way, for the nonbreaking case it is found that α equal to 1.2 and r equal to 5 present better results (both in terms of computational effort and accuracy) for IITM-FNPT and OceanWave3D, respectively. Similarly, for both these two models, where the timestep size is fixed in the time marching scheme, choosing 100 steps per T p improves the CPU time without affecting the accuracy. In an analogous way, it results that the best choices for the time and GMRES tolerance are 10 -4 and 10 -3 , respectively. Finally, comparing the three models one can see how HOS-NWT is able to reach the highest accuracy without requiring significant computational effort. However, OceanWave3D presents competitive CPU time for the range of accuracy it can reach. The results of IITM-FNPT, instead, present both a higher computational effort and lower accuracy. Regarding the breaking case, a first results is that the time marching scheme for OceanWave3D requires a smaller timestep size to obtain satisfactory accuracy. On the contrary, no change in the optimal choices for the time and GMRES tolerance is found with respect to the nonbreaking case. The final comparison between the breaking models suggests that the smoothing technique does not imply an increase in computational effort if compared to the nonbreaking case. The same cannot be said for all the models tested for HOS-NWT. However, the latter can reach significantly greater accuracy. Regarding the three breaking models used for HOS-NWT, the combination of Tian model and hyperviscous filter results in reaching slightly greater accuracy without affecting the computational effort. However, it should be noticed that for all the models the sampling frequency for the error computation is chosen as considerably high (100 steps per T p for the nonbreaking case, 220 steps per T p for the breaking case). Interesting information, especially on the efficiency of the time integration, could be obtained investigating other sampling frequencies.

The present work could be extended in order to better address the efficiency of numerical solvers for ocean wave tanks simulations. Regarding the considered models, the influence of the time integration scheme can be investigated in more details by extending the study of the timestep size (for RK(4)) or time tolerance (for RK4 [START_REF] Bonnefoy | A Fully Spectral 3D Time-Domain Model for Second-Order Simulation of Wavetank Experiments. Part B: Validation; Calibration versus experiments and Sample Application[END_REF]). Interesting conclusions could also be deduced by evaluating the effect of the number of steps after which regridding is applied in IITM-FNPT or of the dealiasing parameter in HOS-NWT. Next, the case of a three-dimensional wavefield could lead to variations in terms of accuracy and computational effort.

Figure 2 :

 2 Figure 2: Vertical nodes configuration

Figure 3 :

 3 Figure 3: Geometric quantities defining ν eddy used in Tian model

Figure 4 :

 4 Figure 4: Example of hyperviscous filter action

Figure 5 :Figure 6 :

 56 Figure 5: LHEEA ocean wave tank

Figure 7 :

 7 Figure 7: Experimental wave elevation (top) and comparison in time domain between the different models over a 20T p zoom (bottom)

Figure 8 :

 8 Figure 8: ISSP error for varying horizontal discretisation and HOS order, time tolerance fixed to 10 -4

Figure 9 :

 9 Figure 9: ϵ 1 error for varying horizontal discretisation and HOS order, time tolerance fixed to 10 -4

Figure 10 :Figure 11 :

 1011 Figure 10: ϵ 2 error for varying horizontal discretisation and HOS order, time tolerance fixed to 10 -4

Figure 16 :Figure 17 :

 1617 Figure 16: Maximum time lag for cross-correlation for varying horizontal discretisation, 15 vertical nodes

Figure 20 :Figure 21 :

 2021 Figure 20: Finite element mesh for α = 1.2 (top) and α = 2 (bottom), 80 horizontal nodes, 9 vertical nodes

Figure 22 :Figure 23 :

 2223 Figure 22: ISSP error for varying mesh resolution, α fixed to 2, 200 steps per T p

Figure 24 :

 24 Figure 24: Maximum time lag for cross-correlation for varying horizontal discretisation, 20 vertical nodes

Figure 25 :

 25 Figure 25: ISSP error for varying timestep size, α fixed to 1.2, 40 nodes per peak wavelength, 15 vertical nodes

Figure 26 :

 26 Figure 26: Experimental wave elevation (top) and comparison in time domain between the different models over a 20T p zoom (bottom)

Figure 27 :

 27 Figure 27: ISSP error for varying horizontal discretisation and HOS order with hyperviscous filter, time tolerance fixed to 10 -4

Figure 31 :

 31 Figure 31: ISSP error for varying time tolerance with Tian hyperviscous model, 40 nodes per peak wavelength and M equal to 5

Figure 32 :

 32 Figure 32: ISSP error for varying mesh resolution with smoothing technique, r fixed to 5, GMRES tolerance 10 -8 , 220 steps per T p

Figure 33 :Figure 34 :

 3334 Figure 33: ϵ 2 error for varying mesh resolution with smoothing technique, r fixed to 5, GMRES tolerance 10 -8 , 220 steps per T p

Figure

Figure 35 :

 35 Figure 35: ISSP error for varying timestep size with smoothing technique, r fixed to 5, 30 nodes per peak wavelength, 12 vertical nodes

Figure 36 :

 36 Figure 36: ISSP error for varying GMRES tolerance with smoothing technique, r fixed to 5, 220 steps per T p , 30 nodes per peak wavelength, 12 vertical nodes

Figure 37 :

 37 Figure 37: CPU time versus ISSP accuracy for OceanWave3D (left) and IITM-FNPT (right) with different values of r and α, respectively

  Figure (39).

Figure 39 :

 39 Figure 39: CPU time versus ISSP accuracy for HOS-NWT (left) and OceanWave3D (right) with different time and GMRES tolerance, respectively

Figure 40 :

 40 Figure 40: CPU time versus ISSP accuracy for the different methods in optimal configuration

Figure 41 :

 41 Figure 41: CPU time versus ISSP accuracy for the different models for all cases tested

Figure 44 :

 44 Figure 44: CPU time versus ISSP accuracy for the different models in optimal configuration

Figure 45 :

 45 Figure 45: CPU time versus ISSP accuracy for the different models for all cases tested

Table 1 and

 1 Table 2 list the parameters that have been presented previously.

		Space discretisation Time discretisation	Specificities
	HOS-NWT	N x	N z	time tolerance	M , dealiasing parameter
	OceanWave3D N x	N z	timestep size	r, GMRES tolerance
	IITM-FNPT	N x	N z	timestep size	α, n. of steps for regridding

Table 1 :

 1 Main parameters for each model

			Model onset	Specificities
	HOS-NWT	Tian Hyperviscous	threshold (0.85) None	β a, k d
	OceanWave3D	Smoothing	γ	None

Table 2 :

 2 Main parameters for each breaking model

Table 3 :

 3 JONSWAP spectra parameters for experimental data

	Case 1	13.6	12.3	1	200	40	17.5	13.95
	Case 2	17	15.5 2.6	250	50	14	11.25

  Figure 14: ISSP error for varying mesh resolution, r fixed to 7, GMRES tolerance 10 -8 , 200 steps per T p

	0.44									0.33
	0.42									0.32
	0.4									0.31
	0.38
										0.3
	0.36
										0.29
	0.34
										0.28
	0.32
										0.27
	0.3
	0.28									0.26
	0.26									0.25
	0.24									0.24
	10	15	20	25	30	35	40	45	50	55	7	8	9	10	11	12	13	14	15
	0.5									0.264
										0.262
	0.45
										0.26
	0.4
										0.258
	0.35									0.256
										0.254
	0.3
										0.252
	0.25
										0.25
	0.2									0.248
	10	15	20	25	30	35	40	45	50	55	9	10	11		12	13	14		15

Figure 15: ISSP error for varying mesh resolution, r fixed to 9, GMRES tolerance 10 -8 , 200 steps per T p

  Figure 28: ϵ 2 error for varying horizontal discretisation and HOS order with hyperviscous filter, time tolerance fixed to 10 -4

	0.15										0.4
	0.14
											0.35
	0.13
											0.3
	0.12
	0.11										0.25
	0.1										0.2
	0.09
											0.15
	0.08
											0.1
	0.07
	0.06										0.05
	15	20		25	30		35	40	45		50	1	2	3	4	5	6	7
	0.31										0.38
	0.3										0.36
	0.29										0.34
	0.28										0.32
	0.27										0.3
	0.26										0.28
	0.25										0.26
	0.24										0.24
	20	22	24	26	28	30	32	34	36	38	40	1	2	3	4	5	6	7
	Figure 29: ISSP error for varying horizontal discretisation and HOS order with Tian model,
	time tolerance fixed to 10 -4

  Figure 30: ISSP error for varying horizontal discretisation and HOS order with Tian hyperviscous model, time tolerance fixed to 10 -4

	0.285							0.38
	0.28
								0.36
	0.275
								0.34
	0.27
	0.265							0.32
	0.26							0.3
	0.255
								0.28
	0.25
								0.26
	0.245
	0.24							0.24
	15	20	25	30	35	40	45	50	1	2	3	4	5	6	7

Table 4 :

 4 Optimal choices for efficiency comparison for HOS-NWT

	report the final configurations

Table 6 :

 6 Optimal choices for efficiency comparison for IITM-FNPT

Table 7 :

 7 . Optimal choices for efficiency comparison for HOS-NWT with hyperviscous filter

	ISSP error Nodes per λ p M
	HOS-NWT; hyperviscous
	0.291	20	4
	0.287	20	5
	0.272	30	5
	0.253	40	4
	0.251	40	5
	0.251	50	5
	ISSP error Nodes per λ p M
	HOS-NWT; Tian
	0.309	20	4
	0.301	20	5
	0.263	30	5
	0.252	40	4
	0.247	40	5

Table 8 :

 8 Optimal choices for efficiency comparison for HOS-NWT with Tian modelAs for the previous case, the best choice of two other parameters (timestep size and GMRES tolerance) should be investigated for OceanWave3D and of the time tolerance for HOS-NWT. First, the timestep size influence for OceanWave3D is presented. As for the

						3
						2.5
						2
						1.5
						1
						0.5
						0
						0.714	0.716	0.718	0.72	0.722	0.724	0.726	0.728	0.73	0.732	0.734
	Figure 42: CPU time versus ISSP accuracy for OceanWave3D with smoothing technique
	with different number of steps per peak waveperiod
	25												3
	20												2.5
	15												2
	10												1.5
	5												1
	0												0.5
	0.71	0.715	0.72	0.725	0.73	0.735	0.74	0.745	0.75	0.755	0.76		0.725	0.726	0.727	0.728	0.729	0.73	0.731	0.732	0.733	0.734
	Figure 43: CPU time versus ISSP accuracy for HOS-NWT with Tian hyperviscous model
	(left) and OceanWave3D with smoothing technique (right) with different time and GMRES
	tolerance, respectively

previous case, the cases studied are those with a length of 0.01s and 0.02s, here approximately corresponding to 110 and 220 steps per T p , respectively. Figure (42) shows the CPU time against the accuracy obtained. Similarly to the nonbreaking case, in this plot and in the following of the same type each point represents one of the choices of discretisation parameters of Tables [START_REF] Ducrozet | A comparative study of two fast nonlinear free-surface water wave models[END_REF] to [START_REF] Gao | Numerical simulation of deterministic freak wave sequences and wavestructure interaction[END_REF]. It is appreciable how a higher computational effort does not always correspond to higher accuracy. This is the direct reflection of the non perfect convergence discussed in Section (5.2.2). Nevertheless, the smaller CPU time required for 110 steps per waveperiod also implies a significant decrease in the accuracy, as could be expected by looking at the ISSP convergence over the timestep size in Figure (35). Therefore, one can conclude that 110 steps per wavepriod are not sufficient to reach a converged numerical result. For this reason, unlike for the nonbreaking case 220 steps per T p will be chosen in the following.

The two last parameters here considered are the time tolerance for HOS-NWT and the GMRES tolerance for OceanWave3D. For the purpose of this analysis, for HOS-NWT only the Tian model with hyperviscous filter is presented, but similar conclusions are found for the other breaking models. Figure (43) shows the efficiency comparison between the different cases tested for HOS-NWT on the left and OceanWave3D on the right.

Similar conclusions to those drawn for the nonbreaking case can be obtained. As it can be noticed, regarding the time tolerance for HOS-NWT, the choice of 10 -4 presents significantly less computational effort. At the same time, the accuracy of the numerical solution dos not decrease greatly; thus this can be considered as an optimal configuration.


Quality Management Systems Alignment in European Pork Chains: A Governance Perspective

Research questions  How to characterize the functions related to the governance of quality?  What objectives do the QMSs fulfill from the quality point of view?  What are the properties of different governance structures? Design Theories  Governance of quality: objectives, functions and levels  Alignment between QMSs: an analytical framework based on a governance perspective Operationalization  Case study approach  In depth interviews Empirical research: application to French case studies  Ham of Bayonne: Protected Geographical Indication applied to dry cured ham in the South West of France  Le Cochon de Bretagne: a collective pig producer's brand in Brittany, West of France  Porcilin/ Saveurs en'OR: pig producer brand coupled with a collective regional brand in the Nord-Pas de Calais region in France  Ham Fleury Michon with omega 3 /Bleu Blanc Coeur: national processor's brand associated with a collective nutritional brand Results  Synthesis based upon case study categories  Background theories about governance of quality  Managerial implications: Global synthesis/Recommendations to public and private stakeholders on cost efficiency, effectiveness of QMSs and control redundancy

Further research

Comparison of the same category of case studies in different countries and country comparison: situation at stake/institutional aspect Rakotonandraina, N., & Sauvée, L. Institut Polytechnique LaSalle Beauvais, France (nalini.rakotonandraina@lasalle-beauvais.fr ; loic.sauvee@lasalle-beauvais.fr)

Abstract: This research proposes an analytical framework in order to assess the alignment of quality management systems (hereafter QMSs) understood as the concrete implementation between different QMSs in European pork chains. The approach uses the concept of governance of quality in the analytical part, defined as a set of functions to solve quality coordination problems. Thus QMSs are analyzed through a governance of quality perspective by considering its fulfilled functions and levels of decision as regard to institutional environment. Four case studies in French pork chains are used to illustrate the analysis. Results contribute in their theoritical part to the institutional economics theory in line with previous works on alignment and at empirical level give managerial perspectives on the alignment of QMSs in terms of efficiency in the one hand (cost, QMSs coherency, quality controls and specific roles of actors) and highligths complementarity and redundancy devices as regard to the QMSs alignment typology in the other hand.

Objectives

 To provide a governance perspective on QMSs in pork chains  To propose a framework to investigate the question of alignment of QMSs within pork chains  To identify the main "alignment situations" between QMSs in real cases This table (4) gives an overview of QMSs alignment typology considering subgroup of QMSs found in the case of Fleury Michon/Bleu Blanc Coeur.

Figure1:Figure3

  Figure1: Operationalization model of governance adaptation Figure2: QMSs classical categorization Table1: QMSs alignment typology principle

Acknowledgement

The authors gratefully acknowledge from the European Community financial participation under the Sixth Framework Programme for Research, Technological Development and Demonstration Activities, for the Integrated Project Q-PORKCHAINS FOOD-CT-2007-036245.

Disclaimer

The views expressed in this publication are the sole responsibility of the author(s) and do not necessarily reflect the views of the European Commission. Neither the European Commission nor any person acting on behalf of the Commission is responsible for the use, which might be made of the information. The information in this document is provided as is and no guarantee or warranty is given that the information is fit for any particular purpose. The user thereof uses the information at its sole risk and liability.