# A comparative study of two fast nonlinear free-surface water wave models

G. Ducrozet, H.B. Bingham, A.P. Engsig-Karup, F. Bonnefoy, P. Ferrant

► To cite this version:

G. Ducrozet, H.B. Bingham, A.P. Engsig-Karup, F. Bonnefoy, P. Ferrant. A comparative study of two fast nonlinear free-surface water wave models. International Journal for Numerical Methods in Fluids, 2012, pp.1818-1834. 10.1002/fld.2672 . hal-01145120

HAL Id: hal-01145120

https://hal.science/hal-01145120

Submitted on 24 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A comparative study of two fast nonlinear free-surface water wave models

Guillaume Ducrozet¹,*,†, Harry B. Bingham², Allan P. Engsig-Karup³, Félicien Bonnefoy¹ and Pierre Ferrant¹

¹Ecole Centrale Nantes - Laboratoire de Mécanique des Fluides, Nantes, France

²Technical University of Denmark, Dept. of Mechanical Engineering, Kgs. Lyngby, Denmark

³Technical University Denmark, Dept. of Informatics and Mathematical Modeling, Lyngby, Denmark

## SUMMARY

This paper presents a comparison in terms of accuracy and efficiency between two fully nonlinear potential flow solvers for the solution of gravity wave propagation. One model is based on the high-order spectral (HOS) method, whereas the second model is the high-order finite difference model OceanWave3D. Although both models solve the nonlinear potential flow problem, they make use of two different approaches. The HOS model uses a modal expansion in the vertical direction to collapse the numerical solution to the two-dimensional horizontal plane. On the other hand, the finite difference model simply directly solves the three-dimensional problem. Both models have been well validated on standard test cases and shown to exhibit attractive convergence properties and an optimal scaling of the computational effort with increasing problem size. These two models are compared for solution of a typical problem: propagation of highly nonlinear periodic waves on a finite constant-depth domain. The HOS model is found to be more efficient than OceanWave3D with a difference dependent on the level of accuracy needed as well as the wave steepness. Also, the higher the order of the finite difference schemes used in OceanWave3D, the closer the results come to the HOS model.

KEY WORDS: hydrodynamics; water waves; high-order finite differences; high-order spectral; Ocean-Wave3D; numerical comparisons

## 1. INTRODUCTION

The study of surface gravity waves (propagation, wave-wave interactions, etc.) is of major interest. The accurate description of wave fields is necessary in the ocean and naval engineering context to determine precisely the nonlinear wave loads acting on structures at sea, for instance.

A wide variety of nonlinear wave models have been developed in the last decades to take care of this problem. Most of them were developed in the framework of potential flow theory considering that ocean wave propagation is essentially irrotational and inviscid (until the point of wave breaking). It has to be noted that attempts have been made to solve the full Euler or Navier-Stokes equations (see Park et al. [1], for instance). However, the computational effort required here restricts the solution to very small scales. At the same time, the efficient and accurate solution of the fully nonlinear potential problem is still very challenging.

This study compares two of the fully nonlinear potential flow models. The first model is a pseudo-spectral model, namely, high-order spectral (HOS), initially developed by West et al. [2]

<footnote>*Correspondence to: G. Ducrozet, Ecole Centrale Nantes - Laboratoire de Mécanique des Fluides ; 1, rue de la Noë - 44321 Nantes, France.</footnote>

<footnote>†E-mail: guillaume.ducrozet@ec-nantes.fr</footnote>
and Dommermuth and Yue [3]. This model has been widely used [4–6], and it has shown its ability to simulate efficiently wave evolution (propagation, nonlinear wave-wave interactions, etc.) in domains of large scale with reasonable computational effort (cf. Ducrozet et al. [7]). The second model, OceanWave3D, is dedicated to the simulation of wave-wave, wave-bottom and wave-structure interactions. It is based on a high-order finite difference solver which has also shown its ability and efficiency (see e.g. Bingham and Zhang [8] (2D) and Engsig-Karup et al. [9] (3D)). Finite differences solution of unsteady free-surface flows is also a widespread numerical method [10–12]. Note that the pseudo-spectral formalism of the HOS model restricts the model to a rectangular domain (in the horizontal plane) and imposes limits on the total variation of the water depth, whereas OceanWave3D is more flexible with regard to the geometry and bathymetry.

Firstly, the two numerical methods which will be compared are presented. The general potential formulation and framework of both methods are reviewed. Then, the numerical approach based on high-order finite differences adopted in OceanWave3D is described. It is followed by the presentation of the HOS numerical scheme. The second and third parts are dedicated to comparing the two models. Firstly, the numerical properties of both methods are presented and particularly the scaling of computational effort and RAM memory use with respect to the space discretization. This leads to a first comparison of the general behaviour of these methods. Then, the two models are compared in terms of efficiency and accuracy. A highly nonlinear wave is computed, and the error of evaluation of the vertical velocity is studied as well as the relative efficiency of both models. The comparison is pursued with the long-time propagation of the previous nonlinear periodic wave. Errors for both models on the phase shift between the computed and analytical solutions are investigated.

## 2. NUMERICAL METHODS

In this section, we present the two models of interest in this comparison. The general framework of the fully nonlinear potential flow solution process is first given. Then, the high-order finite differences model OceanWave3D is detailed. Finally, the HOS model is described. The details of the models are reviewed, whereas the effective comparisons will be presented in Sections 3 and 4. More details on the methods can be found in [7–9].

### 2.1. General framework

We consider a rectangular fluid domain $D$ of horizontal dimensions ($L_x, L_y$). We choose a Cartesian coordinate system with the origin $O$ located at one corner of the domain $D$. The $z$ axis is vertical and oriented upwards, with the level $z = 0$ corresponding to the mean water level. The notation $\mathbf{x}$ stands for the $(x, y)$ vector.

We work in the following with a potential flow formalism. We assume that the fluid is incompressible and inviscid and that the flow is irrotational. With these assumptions (relevant in the context of wave propagation), the velocity $\mathbf{V}$ derives from a potential $V(\mathbf{x}, z, t) = (\nabla, \partial_z) \phi(\mathbf{x}, z, t)$, $\nabla$ representing the horizontal gradient, $\partial_z$ the partial derivative with respect to $z$ and $\phi$ the velocity potential. The continuity equation: $(\nabla, \partial_z) \mathbf{V} = 0$ becomes the Laplace equation in the fluid domain $D$

$$ \Delta\phi = 0 \qquad (1) $$

On the free surface, the elevation is supposed to be single valued and thus described by $z = \eta(\mathbf{x}, t)$. On this boundary, the slip condition and the continuity of the pressure give, respectively, the kinematic and dynamic free-surface conditions. Following Zakharov [13], these fully nonlinear free-surface boundary conditions can be written in terms of surface quantities, namely $\eta$ and the surface potential $\tilde{\phi}(\mathbf{x}, t) = \phi(\mathbf{x}, \eta, t)$

$$ \tilde{\phi}_t = -g\eta - \frac{1}{2}|\nabla\tilde{\phi}|^2 + \frac{1}{2}(1 + |\nabla\eta|^2)\phi_z^2 \qquad (2) $$

$$ \eta_t = (1 + |\nabla\eta|^2) \phi_z - \nabla\tilde{\phi} \cdot \nabla\eta \qquad (3) $$

on $z = \eta(\mathbf{x}, t)$.
To account for the time evolution of the quantities of interest ($\eta$, $\tilde{\phi}$), one only needs to evaluate the vertical velocity at the free surface $\tilde{w}(\mathbf{x}, t) = \phi_z(\mathbf{x}, \eta, t)$. The solution process chosen to evaluate this vertical velocity differs in the two models of interest and is the subject of the next paragraphs. The set of Equations (2) and (3) then provides the time derivatives of the unknowns $\eta$ and $\tilde{\phi}$, which are further used in a time-marching classical fourth-order Runge-Kutta scheme in both methods. In the following, time steps $\Delta t$ are chosen fixing the Courant number $C_r = \frac{C\Delta t}{\Delta x}$ with C standing for the phase velocity of the solved wave and $\Delta x$ the space resolution.

## 2.2. OceanWave3D model

In OceanWave3D, the evaluation of the vertical velocity at the free surface is done by solving the following Laplace problem for the velocity potential $\phi$

$$ \phi = \tilde{\phi}, \quad z = \eta \qquad (4) $$

$$ \nabla^2 \phi + \phi_{zz} = 0, \quad -h \le z < \eta \qquad (5) $$

$$ (\mathbf{n}, n_z).(\nabla, \partial_z)\phi = 0, \quad (\mathbf{x}, z) \in \partial\Omega \qquad (6) $$

with (6) the kinematic boundary conditions on the different boundaries $\partial\Omega$ of the finite domain, including the bottom whose profile can be chosen arbitrarily.

A direct solution of the Laplace problem is computed thanks to the following $\sigma$ transformation $\sigma(\mathbf{x}, z, t) = \frac{z + h(\mathbf{x})}{\eta(\mathbf{x}, t) + h(\mathbf{x})}$. This non-conformal mapping allows the transformation of the physical domain $(\mathbf{x}, z)$ (which evolves during time) onto a fixed computational domain $(\mathbf{x}, \sigma)$.

The Laplace problem expressed in $(\mathbf{x}, \sigma)$-coordinate is then solved by a finite difference method. Arbitrary high-order schemes are used with possibly stretched grid in both directions. These developments are found to be advantageous relative to classical second-order schemes on even grids.

The obtained linear system is solved with a preconditioned iterative method. Note that to impose boundary conditions (lateral walls, bottom), ghost point strategy is used. The original generalized minimal residual (GMRES) method (without restarts) is formally not memory limited, which may lead to excessive memory requirements. Instead, a defect correction method can be employed as an alternative iterative solver at comparable efficiency, and this reduces the memory requirements significantly (see Engsig-Karup [14]).

One initial LU factorization for each preconditioning step is effective in 2D, but for large 3D problems, the direct solution of the preconditioning step is replaced with a multigrid solver (to retain optimal scaling in computational effort and RAM memory use [9]). Note that the multigrid solution of the problem was previously developed by Li and Flemming [12] before the enhancement proposed in Bingham and Zhang [8] or Engsig-Karup et al. [9]. One refers to those two papers for more details about the method and different validations performed.

## 2.3. High-order spectral model

The HOS method was first introduced by West et al. [2] and Dommermuth and Yue [3]. The domain $D$ is assumed periodic in both the $x$ and $y$ horizontal directions, and a slip condition is imposed on the bottom of the domain (if there is one). This allows us to define a spectral basis on which the velocity potential will be expanded. Note that this decomposition on spectral basis enforces the geometry of the computational domain (typically rectangular domains in the horizontal plane with constant depth). The spectral basis for $\phi$ is defined with the basis functions $\psi_{mn}$ which individually satisfy the set of equations: Laplace and periodicity

$$ \phi(\mathbf{x}, z, t) = \sum_{m,n} A_{mn}(t) \psi_{mn}(\mathbf{x}, z) \qquad (7) $$

$$ m_n(\mathbf{x}, z) = \exp(ik_m x) \exp(ik_n y) \frac{\cosh(k_{mn}[z+h])}{\cosh(k_m nh)} \qquad (8) $$
where $k_m = 2m\pi/L_x$, $k_n = 2n\pi/L_y$ and $k_{mn} = \sqrt{k_m^2 + k_n^2}$ are the wavenumbers associated with modes $(m, n) \in \mathbb{Z}^2$. A similar basis based only on the horizontal coordinates is used for the free-surface elevation $\eta$ and surface velocity potential $\tilde{\phi}$. The spectral basis for the potential taken at $z=0$ and the one for surface quantities are adequate for use of Fourier transforms. It has to be noted that a numerical wave tank has also been designed with HOS method, thus accounting for reflective condition in the horizontal directions (see Ducrozet et al. [15, 16]).

The evaluation of the vertical velocity $\tilde{w}(\mathbf{x}, t) = \phi_z(\mathbf{x}, \eta, t)$ in the free-surface boundary conditions is performed thanks to the order-consistent HOS scheme of West et al. [2]. It consists in a double expansion of the potential $\phi$ to solve the Dirichlet problem $\phi(z = \eta) = \tilde{\phi}$. First, the potential $\phi$ is expressed as a truncated power series of components $\phi^{(m)}$ for $m = 0$ to $M$, each component being of magnitude $\eta^m$. Second, the potential taken at the free surface is expanded in a Taylor series about the mean water level $z=0$. Combining these two expansions gives a triangular set of Dirichlet problems for the components that can be solved by means of a spectral method (using efficient fast Fourier transforms). Then by evaluating the vertical derivative of the potential and using the same kind of double expansion, we obtain the vertical velocity on the free surface from the components $\phi^{(m)}$.

The products involving $\nabla\eta$ and $W$ in free-surface boundary conditions are evaluated thanks to the order-consistent formulation of [2]. Note that in this process, it is crucial to perform a careful dealiasing to preserve the method convergence and accuracy for waves close to Stokes’ limit [17].

This HOS model demonstrates its accuracy and efficiency (see [7, 17]). It has been validated on several test cases: propagation of highly nonlinear regular and irregular waves and study of wave-wave interactions (with application to the formation of freak waves).

### 3. SCALING OF COMPUTATIONAL EFFORT AND MEMORY USE

In this section is provided a preliminary study of the models of interest. The scaling of computational effort as well as the RAM memory use is studied with respect to the number of points in the simulation with both models (firstly with OceanWave3D code and then with the HOS model).

#### 3.1. OceanWave3D

Figure 1 gives the scaling in terms of RAM memory use and computational effort for the Ocean-Wave3D code with MultiGrid preconditioning for typical 3D computations (nonlinear periodic wave studied in Section 4.1) run on an Intel Xeon E5520 2.27 GHz. The total RAM memory used (left part of the figure) is given as well as the CPU time per defect correction (DC) iteration (right part

<figure><img src="image_2.png" /><figcaption>Figure 1. Scaling RAM memory use and computational effort—3D simulation, sixth order.</figcaption></figure>
of the figure) as a function of the number of points in the computational domain $N = N_x N_y N_z$. The dashed lines represent the expected linear scaling. Note that a typical sixth-order finite difference schemes is used in this study, but results are similar to any order of the finite differences (FD) scheme.

It appears from this figure that the memory requirement evolves linearly with the number of points in the domain. At the same time, the computational effort scales are also perfectly linear (even for quite a large number of points in the domain: here up to $2.4 \times 10^6$ points). This study is made on 3D configuration because the *OceanWave3D* model has been especially designed for this and this is the most demanding test case (to achieve this linear scaling). This is the multigrid preconditioning strategy which allows to retain the linear scaling of both computational effort and memory requirement (see Engsig-Karup *et al.* [9]). Note that these linear scalings are also obtained for 2D computations.

## 3.2. High-order spectral model

Figure 2 presents the same plots as in the previous paragraph with the scaling of RAM memory use (left part) and computational effort per Runge-Kutta step (right part) with respect to the number of points ($N_x N_y$) for the HOS method. This figure was obtained for a typical 3D simulation with a HOS order fixed to $M = 5$. The right part of this figure is slightly different from the one obtained with *OceanWave3D*, which plots the computational effort per DC iteration. Furthermore, note that the HOS model takes full advantage of the Laplace equation and thus has unknowns only on the free surface: the number of points involved in this figure is $N = N_x N_y$ (the vertical direction is not discretized). However, the global behaviour of each of the methods may be compared, keeping these features in mind.

Results indicate that the memory requirements scale linearly with the number of points in the domain. Dealing with the computational effort, it appears that the HOS model scales with $N \log(N)$, as expected. The $N$ scaling is included on this figure to point out that it is also close to the behaviour of the HOS model. The same behaviour is observed for typical 2D simulations.

## 3.3. Discussion

In order to compare both models in terms of efficiency and accuracy from the previous results, one has to recall the details of each model:

* Firstly, as noted before, the number of points taken into account in the two models is different because the HOS method solves the problem on the free surface only, whereas in the *OceanWave3D* code, the vertical direction is also discretized. The number of points $N$ in Figures 1

<figure><img src="image_3.png" /><figcaption>Figure 2. Scaling RAM memory use and computational effort—3D simulations, $M = 5$.</figcaption></figure>
and 2 is thus not equivalent. However, it has been shown [8, 9] that the high-order finite dif-
ferences scheme of OceanWave3D requires few points in the vertical direction to achieve good
accuracy.

* Secondly, the results will be strongly dependent on the choice of the order of the finite difference scheme and the DC or GMRES tolerance for OceanWave3D model and dependent on the order of nonlinearity *M* for the HOS model. This will influence both RAM memory use and computational effort. However, note that the *O(N)* and *O(N log(N))* scalings are conserved with all choices of parameters in both methods.

* Finally, in Figures 1 and 2 is given the evolution of computational effort respectively per DC iteration and per Runge-Kutta step. The comparison is thus subtle even if the solution with OceanWave3D is obtained in $O(10)$ DC or GMRES iteration per Runge-Kutta step almost independent of number of points in the domain and physics of the studied phenomenon.

Consequently, memory use appears to be equivalent between both models, whereas computa-
tional effort probably admits a break-even point between both models (*O*(N) versus *O*(N log(N))
behaviour). However, previously listed details indicate that it will be complicated to esti-
mate/determine this break-even point for 2D or 3D problems (number of points different in both
approaches and results dependent on the numerical parameters).

Next, Section 4 proposes a comparison of the efficiency and accuracy of both models. It is based
on a study of numerical methods which have to reach a given level of accuracy for the solution of a
standard 2D test case.

4. EFFICIENCY AND ACCURACY

The efficiency and accuracy of both models is now analysed. The approach chosen is to study OceanWave3D and HOS models on a typical test case (propagation of nonlinear periodic wave). A given level of accuracy has to be achieved, and the corresponding computational effort needed with both models is studied. A first test case is studied in which the vertical velocity at the free surface of a 2D nonlinear periodic wave is computed at one time step and compared with the stream function reference solution. A second test case is then studied where the same 2D nonlinear periodic wave is now propagated (i.e. advanced in time) with the two models, and phase shift observed after a long propagation time (1000 periods) is monitored.

4.1. *First test case—vertical velocity of a 2D periodic wave*

Recall that both methods consider the fully nonlinear free-surface boundary conditions formulated,
following Zakharov [13], Equations (3) and (2). In order to advance in time the unknowns on the
free surface $\eta$ and $\tilde{\phi}$, one needs to evaluate the vertical velocity $\tilde{w}$. This first test case proposes to
study, for a nonlinear periodic wave, the error made on this evaluation of the vertical velocity at a
given time step. We therefore consider for now only errors due to the spatial solution (time-stepping
errors being considered in Section 4.2).

This wave is described by its height $H$ (or equivalent amplitude $a = H/2$), wavelength $\lambda$ (or
wavenumber $k = 2\pi/\lambda$) and water depth $h$. The nonlinear wave we choose here has the following
characteristics: $H/\lambda = 0.0955 = 67.5\%$ of limiting steepness, that is $ka = 0.3$ and $kh = 2\pi$ (i.e.
deep water). The computational domain is periodic in both models and of length $\lambda$. The initial con-
dition (free-surface elevation $\eta$ and velocity potential $\tilde{\phi}$) is known from stream function theory [18].
The vertical velocity obtained using those $\eta$ and $\tilde{\phi}$ with both models is compared with the stream
function theory assumed to be exact (convergence up to machine accuracy achieved).

The results obtained with *OceanWave3D* are first presented. Then, the computations are per-
formed with the HOS method. Finally, a comparison between both methods is proposed.

4.1.1. OceanWave3D. In this study, regular grid is used in the horizontal direction with periodic boundary conditions. In the vertical direction, a clustering of the nodes towards the free surface is applied to enhance the efficiency of the solution (similar to the one used in Bingham and Zhang [8], i.e. a cosine spaced grid). The tolerance of the GMRES method is fixed to $10^{-15}$ to ensure a correct
convergence with respect to the number of points. Note that for the effective comparisons, it will be adjusted to the level of accuracy needed (see Section 4.1.3).

Figures 3, 4 and 5 present the influence of the discretization $N_x$ and $N_z$ on the error of evaluation of the vertical velocity $\tilde{w}$. The error is defined as the maximum of $|\tilde{w} - \tilde{w}_{ref}|$, with $\tilde{w}_{ref}$ standing for the reference solution obtained from stream function theory [18]. In these figures, different orders of the finite difference scheme are tested (respectively, fourth, sixth and eighth orders). The left part of the figure gives the convergence with respect to the vertical discretization $N_z$, whereas the right part deals with the horizontal discretization $N_x$. On these plots is given in dashed line the theoretical convergence rate of each one of the finite differences schemes used (respectively, fourth, sixth or eighth order).

These plots indicate the effective convergence of the solution with respect to the number of points in the domain (in both directions $N_x$ and $N_z$). An asymptotic slope is obtained with the finest discretizations ($N_x$ = 1024 for left plots and $N_z$ = 129 for right plots). This slope is in perfect agreement with the expected convergence rate.

It is also noticeable that the level of accuracy achieved with the models at a given number of points is greatly improved with high-order finite differences schemes. In addition, it must be emphasized

<figure><img src="image_4.png" /><figcaption>Figure 3. Evolution of error on $\tilde{w}$ with respect to $N_x$ and $N_z$ in OceanWave3D: GMRES tol = $10^{-15}$, fourth-order FD scheme.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Evolution of error on $\tilde{w}$ with respect to $N_x$ and $N_z$ in OceanWave3D: GMRES tol = $10^{-15}$, sixth-order FD scheme.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5. Evolution of error on $\tilde{w}$ with respect to $N_x$ and $N_z$ in OceanWave3D: GMRES tol = $10^{-15}$, eighth-order FD scheme.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Evolution of error on $\tilde{w}$ with respect to finite difference order in OceanWave3D: GMRES tol = $10^{-15}$, $N_z = 65$ (left) and $N_x = 512$ (right).</figcaption></figure>

that the computational effort required to go from fourth-order to eighth-order scheme at a given number of points is roughly an increase of 25%. Then, using high-order finite differences schemes possibly leads to a great improvement of computational effort. Thus, two types of convergence may be used with OceanWave3D code, namely *h-type* and *p-type*, corresponding to convergence relative to the discretization and to the order of finite differences scheme, respectively (which order is left arbitrary).

Figure 6 presents the convergence plot with respect to the order of the finite difference scheme. It appears, as expected, that one achieves an exponential convergence rate (note the log scale in y direction versus the linear scale in x). The left part of the figure depicts the behaviour at different discretizations in horizontal plane $N_x$ with $N_z = 65$. The right part gives the results for different discretizations in vertical plane $N_z$ with $N_x = 512$.

4.1.2. *High-order spectral model.* The same convergence study done previously with Ocean-Wave3D is performed with the HOS model. Figure 7 shows the evolution of the error on the vertical velocity $\tilde{w}$ with respect to the number of points/modes $N_x$ and the HOS order $M$. These are the
two parameters involved in the convergence process of the HOS method. Because one deals with pseudo-spectral method, an exponential convergence rate is expected with respect to the number of points.

The proposed 3D-schematization enables a global view of both convergence with respect to number of points *N* and HOS order *M*. This is a typical plot for such studies obtained with HOS method. It has to be noted that the error is given in log scale (vertical axis), whereas *N* and *M* are plotted in a linear scale in this figure. Two planar surfaces are sketched on Figure 7 (right), which confirms that (i) at high HOS order *M*, the convergence rate is exponential with respect to the number of points (surface A) and (ii) at high number of points *N*, the convergence rate is also exponential with respect to the HOS order *M* (surface B).

Figure 8 is a 2D representation of the 3D plots in Figure 7. The convergence with respect to the number of points on the free surface *N* is given at different HOS orders *M* (left part) as well as the one with HOS order *M* at different numbers of points *N*. This allows a clearer comparison with OceanWave3D results as well as more accurate view of the level of accuracy reached with the HOS model.

The exponential convergence rate appears clearly (straight lines obtained in semi-log plot). This is deduced as a function of number of points (left part of the plot) and as a function of the HOS order

<figure><img src="image_8.png" /><figcaption>Figure 7. Evolution of error on $\tilde{w}$ with respect to $N_x$ and $M$ in HOS model.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Evolution of error on $\tilde{w}$ with respect to $N_x$ and $M$ in HOS model.</figcaption></figure>
(right part of the plot). One obtains the asymptotic behaviours of the HOS model for this nonlinear periodic wave test case:

$$ \text{Error on } \tilde{w} = 10^{(-0.161N - 1.15)} \quad \text{for } M > 0.3N $$

$$ \text{Error on } \tilde{w} = 10^{(-0.508M - 0.718)} \quad \text{for } N > 3.3M $$

### 4.1.3. Comparison of models.
Results with both methods have been presented in Sections 4.1.2 and 4.1.1 and will be now compared. The comparisons are realized with the following procedure to ensure that each model is used in an optimal configuration (at least the closest to this optimal one):

* From the convergence plots of OceanWave3D (Figures 3, 4 and 5), one chooses an optimal couple $N_x/N_z$ for each finite differences scheme.

* This optimum is expected to be the intersection between horizontal flat parts (at given $N_x$ or $N_z$) and the part which follows the asymptotic convergence rate (−4, −6 or −8 slope, respectively).

* Look at both plots Error = $f(N_z)$ and Error = $f(N_x)$ to determine the optimal choice for the couple $N_x/N_z$.

* We note that Figures 3, 4 and 5 are obtained with a tolerance of the GMRES method fixed to $10^{-15}$. However, this can be adjusted to the level of accuracy that has to be reached. For each couple $N_x/N_z$, one determines the largest tolerance which gives the same accuracy as the one obtained with tolerance $10^{-15}$.

* From the convergence plots of the HOS method (Figures 7 and 8), one chooses an optimal couple $N_x/M$ with a similar procedure as the one used for OceanWave3D.

* The optimum is expected to be the intersection of the planar surfaces observed in Figure 7.

* At each level of accuracy needed, a couple $N_x/M$ corresponds along this intersection.

Table I presents for different accuracies the optimal choices of parameters to make the comparison between the models ($N_x$, $N_z$ and GMRES tolerance for OceanWave3D and $N_x$ and $M$ for the HOS method).

Figure 9 presents the comparison between the HOS method and OceanWave3D with the set of parameters described in Table I. The total CPU time per period of propagation $T$ is presented as a function of the accuracy on the vertical velocity $\tilde{w}$. The time integration is the same for both models (i.e. Runge-Kutta 4 method with same Courant number). Different orders of the finite differences scheme in OceanWave3D are presented and compared with the HOS method. One can see that the HOS model is always faster than OceanWave3D. As a result, we chose to plot the following results in terms of the efficiency ratio of OceanWave3D versus HOS. This is simply defined as the ratio, at given accuracy, between the computational effort required by the OceanWave3D computation and that by the HOS. Figure 9 (right) shows this efficiency ratio for the same results as in Figure 9 (left).

One may note on the left part of the figure that the HOS model is more efficient than OceanWave3D for the solution of the propagation of nonlinear periodic waves over flat bottom. At a given order of the finite differences scheme, the difference between HOS and OceanWave3D is greater as the level of accuracy required on $\tilde{w}$ is improved. This is not surprising considering the exponential convergence rate of the pseudo-spectral method. By the way, interesting features can be deduced from the previous figure. For instance, comparing sixth-order and eighth-order schemes informs us about the general behaviour of the numerical model. Break-even points (in terms of efficiency) exist (as the one for accuracy $\approx 6.10^{-6}$) which indicate that typically, the lower orders finite differences schemes will be more efficient for solution with low level of accuracy needed. Other break-even points may be observed between fourth and sixth or eighth order for very low accuracy.

At the same time, the higher the order is, the closer to the pseudo-spectral method result it will be at high accuracy levels. Note that this feature has been noted by Fornberg [19], which compares in a general formalism pseudo-spectral methods and finite differences schemes. His conclusion is that for periodic cases, the two approaches are equivalent, comparing a pseudo-spectral method and a finite difference solution involving all points in the domain $N_x$ (consequently of order $N_x - 1$).
<table><caption>Table I. Optimal choices for comparison of efficiency (periodic waves ka = 0.3 and kh = 2π).</caption><thead><tr><th>Accuracy on <i>w&#771;</i></th><th>N<sub>x</sub></th><th>N<sub>z</sub></th><th>Tol.</th></tr></thead><tbody><tr><td colspan="4"><i>OceanWave3D</i> 4th order</td></tr><tr><td>2.3 &times; 10<sup>-3</sup></td><td>16</td><td>7</td><td>10<sup>-7</sup></td></tr><tr><td>1.6 &times; 10<sup>-3</sup></td><td>32</td><td>7</td><td>10<sup>-8</sup></td></tr><tr><td>1.4 &times; 10<sup>-4</sup></td><td>32</td><td>17</td><td>10<sup>-8</sup></td></tr><tr><td>9.4 &times; 10<sup>-6</sup></td><td>64</td><td>33</td><td>10<sup>-11</sup></td></tr><tr><td>6.1 &times; 10<sup>-7</sup></td><td>128</td><td>65</td><td>10<sup>-12</sup></td></tr><tr><td>3.9 &times; 10<sup>-8</sup></td><td>256</td><td>129</td><td>10<sup>-13</sup></td></tr><tr><td colspan="4"><i>OceanWave3D</i> 8th order</td></tr><tr><td>1.9 &times; 10<sup>-3</sup></td><td>10</td><td>9</td><td>10<sup>-7</sup></td></tr><tr><td>4.3 &times; 10<sup>-4</sup></td><td>16</td><td>9</td><td>10<sup>-9</sup></td></tr><tr><td>3.5 &times; 10<sup>-5</sup></td><td>32</td><td>9</td><td>10<sup>-9</sup></td></tr><tr><td>2.1 &times; 10<sup>-7</sup></td><td>64</td><td>17</td><td>10<sup>-11</sup></td></tr><tr><td>1.3 &times; 10<sup>-9</sup></td><td>128</td><td>33</td><td>10<sup>-13</sup></td></tr><tr><td>6.0 &times; 10<sup>-12</sup></td><td>256</td><td>65</td><td>10<sup>-15</sup></td></tr><tr><td colspan="4"><i>OceanWave3D</i> 6th order</td></tr><tr><td>1.0 &times; 10<sup>-3</sup></td><td>16</td><td>7</td><td>10<sup>-6</sup></td></tr><tr><td>6.0 &times; 10<sup>-5</sup></td><td>32</td><td>9</td><td>10<sup>-8</sup></td></tr><tr><td>2.5 &times; 10<sup>-6</sup></td><td>64</td><td>17</td><td>10<sup>-10</sup></td></tr><tr><td>4.7 &times; 10<sup>-8</sup></td><td>128</td><td>33</td><td>10<sup>-12</sup></td></tr><tr><td>7.7 &times; 10<sup>-10</sup></td><td>256</td><td>65</td><td>10<sup>-13</sup></td></tr><tr><td>1.3 &times; 10<sup>-11</sup></td><td>512</td><td>129</td><td>10<sup>-15</sup></td></tr></tbody></table><table><thead><tr><th>Accuracy on <i>w&#771;</i></th><th>N<sub>x</sub></th><th>M</th></tr></thead><tbody><tr><td colspan="3"><i>HOS</i></td></tr><tr><td>3.4 &times; 10<sup>-3</sup></td><td>8</td><td>5</td></tr><tr><td>1.9 &times; 10<sup>-4</sup></td><td>16</td><td>7</td></tr><tr><td>5.9 &times; 10<sup>-6</sup></td><td>27</td><td>9</td></tr><tr><td>8.0 &times; 10<sup>-7</sup></td><td>32</td><td>11</td></tr><tr><td>1.6 &times; 10<sup>-9</sup></td><td>48</td><td>17</td></tr><tr><td>3.5 &times; 10<sup>-12</sup></td><td>64</td><td>23</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>Figure 9. Comparison of CPU time with respect to accuracy on vertical velocity $\tilde{w}$ for OceanWave3D and HOS method.</figcaption></figure>
The previous tests have been performed with the study of a nonlinear periodic wave over flat bottom with the following parameters: $H/λ = 0.0955 = 67.5\%$ of limiting steepness, that is $ka = 0.3$ and $kh = 2π$ (i.e. deep water). It seems interesting to ensure that the previous efficiency comparison is representative of a large range of $kh$ and $ka$ values, that is that there is no dramatic change for some regimes.

*OceanWave3D* and HOS models have been shown to be able to treat efficiently the problem of wave propagation over a large range of wave steepness $ka$ and relative water depth $kh$. In the following, one studies the efficiency ratio introduced previously as a function of these two parameters $ka$ and $kh$. Note that we fix the order of the finite differences scheme used in *OceanWave3D* to 8 for this study (because this is assumed not to play a key role in the study of interest here).

#### 4.1.4. Influence of ka.
First, we fix the relative water depth to $kh = 2π$ (i.e. deep water case), and we look at the evolution of the efficiency ratio between eighth-order *OceanWave3D* and the HOS model for different wave steepnesses $ka = 0.1, ka = 0.3$ (previous case) and $ka = 0.4$ (i.e. $H/λ = 0.0127 = 91\%$ of limiting steepness). Results are presented in Figure 10. Then, it seems clear that when one increases the steepness, the relative efficiency of *OceanWave3D* is better (i.e. closer to HOS) up to a ratio 6–7 for $ka = 0.4$, whereas for $ka = 0.1$ the ratio is at least 50. However, the same behaviour is observed at the different wave steepnesses: the difference between HOS and *OceanWave3D* is greater as the accuracy required on $\tilde{w}$ is improved. Thus, only the mean level of this efficiency ratio is changing with steepness of the studied wave. One may note that for the steepest case ($ka = 0.4$) which is at 91% of limiting steepness, HOS model is close to its stability limit.

#### 4.1.5. Influence of kh.
The same study is provided with the analysis of the influence of the relative water depth, $kh$. In this concern, the ratio to steepness limit is fixed to $H/λ \simeq 67.5\%$. The different test cases are consequently as follows: (i) $kh = 2π$, $H/λ = 0.0955$; (ii) $kh = 1$, $H/λ = 0.0732$; and (iii) $kh = 0.5$, $H/λ = 0.04$. Corresponding results are presented in Figure 11. In contrast to the $ka$ parameter, it appears here that the relative efficiency between HOS and *OceanWave3D* is not influenced by the relative depth of the simulated wave. Both models are able to simulate wave in relatively shallow water with high accuracy.

To conclude on these first results, a first comparison of the efficiency of *OceanWave3D* code and the HOS model has been performed. The error of evaluation of the vertical velocity of a nonlinear periodic wave over a flat bottom has been studied. One looks for optimal choices of numerical

<figure><img src="image_11.png" /><figcaption>Figure 10. Influence of $ka$ on efficiency ratio between *OceanWave3D* and HOS.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11. Influence of $kh$ on efficiency ratio between OceanWave3D and HOS.</figcaption></figure>

parameters in order to reach a given level of accuracy for both methods. Then, the HOS method appears to be more efficient than OceanWave3D code with an efficiency ratio which appears to be highly dependent on the order of accuracy one wants to achieve (HOS model is more efficient as the accuracy required on $\tilde{w}$ is improved). Furthermore, it appears that the higher the order of the finite differences scheme in OceanWave3D is, the closer to the pseudo-spectral method result it will be at high accuracy levels. Those results are dependent on the computed wave and particularly its steepness (relative water depth appears to have no influence): the steeper the wave is, the closer are efficiencies between models.

## 4.2. Second test case—long-time propagation and phase shift

The nonlinear periodic wave used in the first part of Section 4.1 is studied ($H/\lambda = 0.0955 = 67.5\%$ of limiting steepness, that is $ka = 0.3$ and $kh = 2\pi$). After the accuracy on the vertical velocity at one time step, one studies here the phase shift observed after a long propagation time (typically 1000 periods of propagation) compared with results from stream function theory [18]. This allows now to consider errors due to the spatial solution as well as time stepping. Such a study was proposed by Fructus et al. [20]. Particularly, the stream function theory gives the initial condition as well as the nonlinear period, which is different from the linear period $T_{\text{lin}} = \frac{2\pi}{\omega_{\text{lin}}} = \frac{2\pi}{\sqrt{gk \tanh(kh)}}$.

The results achieved with OceanWave3D code are presented in a first section, followed by the ones with the HOS method. Finally, the two numerical approaches are compared.

### 4.2.1. OceanWave3D
Following the results obtained in subsection 4.1.3, the optimal configurations for this study of phase shift are assumed to be the same as the ones deduced during the analysis of the vertical velocity $\tilde{w}$. Consequently, the most efficient set-up of numerical parameters in terms of computational effort are deduced from Figure 9 for the OceanWave3D code. Table II presents the results of the phase shift observed after 1000 periods of propagation as well as the relative error on the amplitude of the wave $\epsilon_a$. The latter is evaluated as $\epsilon_a = \frac{|\eta_{\text{max}}(t = 1000T) - \eta_{\text{max}}(t = 0)|}{\eta_{\text{max}}(t = 0)}$. The Courant number $C_r$ is chosen equal to 0.64.

The results presented in Table II indicate that the behaviour with respect to accuracy is conserved, and it is expected that the parameters chosen are the optimal ones for this long propagation study.

At the same time, the convergence properties of OceanWave3D code have been checked with respect to the different parameters involved in the solution: $N_x$, $N_z$, order of the finite differences
<table><caption>Table II. OceanWave3D: Error of phase shift and wave amplitude $\epsilon_a$ after 1000 periods of propagation, $C_r = 0.64$ (periodic waves $ka = 0.3$ and $kh = 2\pi$).</caption><thead><tr><th>N<sub>x</sub></th><th>N<sub>z</sub></th><th>FD order</th><th>Tol.</th><th>CPU Time (s)</th><th>Phase Shift (deg)</th><th>&epsilon;<sub>a</sub></th></tr></thead><tbody><tr><td>16</td><td>7</td><td>6</td><td>10<sup>-6</sup></td><td>2.3 &times; 10<sup>1</sup></td><td>2.9 &times; 10<sup>3</sup></td><td>2.0 &times; 10<sup>-1</sup></td></tr><tr><td>32</td><td>9</td><td>6</td><td>10<sup>-8</sup></td><td>1.4 &times; 10<sup>2</sup></td><td>2.3 &times; 10<sup>2</sup></td><td>3.3 &times; 10<sup>-2</sup></td></tr><tr><td>32</td><td>9</td><td>8</td><td>10<sup>-9</sup></td><td>2.8 &times; 10<sup>2</sup></td><td>7.8 &times; 10<sup>1</sup></td><td>1.5 &times; 10<sup>-2</sup></td></tr><tr><td>64</td><td>17</td><td>8</td><td>10<sup>-11</sup></td><td>1.6 &times; 10<sup>3</sup></td><td>3.8</td><td>1.0 &times; 10<sup>-3</sup></td></tr><tr><td>128</td><td>33</td><td>8</td><td>10<sup>-13</sup></td><td>2.1 &times; 10<sup>4</sup></td><td>4.4 &times; 10<sup>-1</sup></td><td>1.3 &times; 10<sup>-4</sup></td></tr><tr><td>256</td><td>65</td><td>8</td><td>10<sup>-15</sup></td><td>2.8 &times; 10<sup>5</sup></td><td>9.7 &times; 10<sup>-2</sup></td><td>2.5 &times; 10<sup>-5</sup></td></tr></tbody></table>

scheme, GMRES tolerance (Tol) and Courant number ($C_r$). Figure 12 presents results of the influence of these several parameters on the phase shift after 1000 periods of propagation. All parameters are fixed except the one of interest on each of the subfigures. OceanWave3D thus appears to properly converge with respect to all numerical parameters involved in such a computation.

### 4.2.2. High-order spectral model.
The same study as in the previous paragraph is made with the HOS model, and corresponding optimal results are presented in Table III.

One may note that similar levels of accuracy are obtained with both methods, which confirms that the optimal choices of parameters are probably conserved between $\tilde{w}$ and phase-shift analyses.

We indicate that a similar study has been made in Bonnefoy et al. [17] with slightly better results in terms of phase shift. This is due to the time integration process which is chosen to be fourth-order Runge-Kutta scheme here for consistency with OceanWave3D, instead of fourth-order Runge-Kutta Cash-Karp scheme with adaptive step size in [17]. Furthermore, some choice of numerical parameters are also slightly different in [17].

<figure><img src="image_13.png" /><figcaption>Figure 12. Convergence checking with respect to $N_x$, $N_z$, order of FD scheme, tolerance of GMRES and Courant number $C_r$ for OceanWave3D.</figcaption></figure>
<table><caption>Table III. HOS: Error of phase shift and wave amplitude $\epsilon_a$ after 1000 periods of propagation, $C_r = 0.64$ (periodic waves $ka = 0.3$ and $kh = 2\pi$).</caption><thead><tr><th>N<sub>x</sub></th><th>M</th><th>CPU Time (s)</th><th>Phase Shift (deg)</th><th>&epsilon;<sub>a</sub></th></tr></thead><tbody><tr><td>8</td><td>5</td><td>2.5</td><td>4.1 &times; 10<sup>2</sup></td><td>1.1 &times; 10<sup>-2</sup></td></tr><tr><td>16</td><td>7</td><td>9.2</td><td>3.5 &times; 10<sup>1</sup></td><td>2.0 &times; 10<sup>-3</sup></td></tr><tr><td>27</td><td>9</td><td>4.5 &times; 10<sup>1</sup></td><td>3.6</td><td>1.9 &times; 10<sup>-4</sup></td></tr><tr><td>32</td><td>11</td><td>8.6 &times; 10<sup>1</sup></td><td>1.5</td><td>8.5 &times; 10<sup>-5</sup></td></tr><tr><td>48</td><td>17</td><td>6.7 &times; 10<sup>2</sup></td><td>1.0 &times; 10<sup>-1</sup></td><td>7.0 &times; 10<sup>-6</sup></td></tr><tr><td>64</td><td>23</td><td>2.5 &times; 10<sup>3</sup></td><td>5.3 &times; 10<sup>-2</sup></td><td>3.5 &times; 10<sup>-6</sup></td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>Figure 13. Comparison of CPU time with respect to accuracy on phase shift for OceanWave3D and HOS method.</figcaption></figure>

### 4.2.3. Comparison of models.
As a final comparison, the relative efficiency between OceanWave3D code and HOS model is presented in Figure 13.

In this plot, the line with symbols gives the efficiency of OceanWave3D code with optimal parameters defined from the study with vertical velocity $\tilde{w}$ compared with HOS model (Tables II and III). The symbols give the relative efficiency of other tests performed during the convergence study (i.e. different choices of parameters). This indicates that one is fairly close to the optimal choices in terms of parameters because the curve is approximatively the lower limit of symbol cloud (except the point at best accuracy $\approx 8.8 \times 10^{-2}$).

The same behaviour as previously observed during the comparison on the vertical velocity $\tilde{w}$ is obtained. When the required accuracy is refined, the HOS method is more efficient. It appears that, comparing with Figure 9, the relative efficiency of OceanWave3D with respect to HOS is greatly reduced (a factor around 10 when dealing with this long-time propagation).

## 5. CONCLUSION

An extensive comparison between two models is presented. Firstly, both methods solving the general problem of nonlinear potential flow problems are presented, namely, OceanWave3D, based on high-order finite differences schemes, and the pseudo-spectral HOS method.

A first analysis of both models is provided with respect to the scaling in terms of RAM memory use and computational effort. Several conclusions can be given: the RAM memory use evolves like the number of points $N$, whereas the computational effort evolves like $N$ and $N \log(N)$, respectively, for the OceanWave3D and the HOS model. Memory use appears to be of the same order of magnitude between both models, whereas computational effort probably admits a break-even point ($O(N)$ versus $O(N \log(N))$ behaviour). However, it will be complicated to estimate/determine this
break-even point for 2D or 3D problems in terms of computational effort as a function of the num-
ber of points. It appears indeed that the comparison is not straightforward because of the following:
(i) the number of points taken into account for *OceanWave3D* and in the HOS method is different
because the HOS method is discretized on the surface only, whereas in the *OceanWave3D* code the
vertical direction is also discretized, and (ii) the results will be strongly dependent on the choice of
the order of the finite difference scheme and the GMRES tolerance for *OceanWave3D* model and of
the order of nonlinearity *M* for the HOS model.

Afterwards, the simple case of the solution of the propagation of a nonlinear periodic wave over a
flat bottom is studied. A first comparison of the efficiency of *OceanWave3D* and the HOS model has
been performed regarding the error on the evaluation of the vertical velocity. We selected optimal
choices of numerical parameters in order to reach a given level of accuracy for both methods. Then,
the HOS method appears to be more efficient than *OceanWave3D* code with an efficiency ratio which
is highly dependent on the order of accuracy one wants to achieve. Break-even points (in terms of
efficiency) exist, which indicate that typically, the lower-order finite differences schemes will be
more efficient for solution with low level of accuracy needed. At the same time, the higher the order
is, the closer to the pseudo-spectral method results it will be at high accuracy levels. Furthermore,
it appears that these results are independent of the relative depth of the periodic wave simulated.
However, the relative efficiency between both models is dependent on the wave steepness of the
computed wave (the steeper the wave is, the closer are the efficiencies between both models).

The study of the phase shift observed after a long propagation time (typically 1000 wave periods) compared with the results from stream function theory [18] confirms the previous conclusions. It appears that on this long-time propagation case, the relative efficiency of *OceanWave3D* code compared with HOS model is reduced by a factor around 10 compared with the previous study (i.e. efficiency ratio at least 40 in the range of accuracy studied).

The proposed benchmarks in this paper are performed with 2D waves. However, it is assumed that results would be equivalent when dealing with 3D sea states. Indeed, both methods retain their numerical features (accuracy and efficiency) for 3D configurations. Consequently, differences between models may be slightly modified, but the main conclusions should be the same.

However, one may keep in mind that the flexibility of the *OceanWave3D* model in terms of the
geometry of the computational domain is not exploited here. Indeed, it has been designed to treat the
problem of wave–wave, wave–bottom and wave–structure interaction. Nevertheless, it is interesting
to point out its relative efficiency compared with a method designed to solve the wave propaga-
tion only in rectangular domain. Furthermore, it has to be noted that extension of HOS to variable
bathymetry exists in the literature. An update of this comparison for configurations with wave–wave
and wave–bottom interactions would be interesting.

ACKNOWLEDGEMENTS

The authors wish to thank Direction Générale de l'Armement (project PDE 08C0039) and the Danish
Research Council for Technology and Production Sciences (grant no 09-070032) for partially funding this
work, and the Danish Center for Scientific Computing for supercomputing resources.

REFERENCES

1. Park JC, Kim MH, Miyata H, Chun HH. Fully nonlinear numerical wave tank (NWT) simulations and wave run-up prediction around 3-D structures. *Ocean Engineering* 2003; **30**(15):1969–1996.

2. West B, Brueckner K, Janda R, Milder M, Milton R. A new numerical method for surface hydrodynamics. *Journal of Geophysical Research* 1987; **92**:11 803–11 824.

3. Dommermuth D, Yue D. A high-order spectral method for the study of nonlinear gravity waves. *Journal of Fluid Mechanics* 1987; **184**:267–288.

4. Tanaka M. Verification of Hasselmann’s energy transfer among surface gravity waves by direct numerical simulations of primitive equations. *Journal of Fluid Mechanics* 2001; **444**:199–221.

5. Toffoli A, Onorato M, Bitner-Gregersen E, Monbaliu J. Development of a bimodal structure in ocean wave spectra. *Journal of Geophysical Research* 2010; **115**:C03006. DOI: 10.1029/2009JC005495.

6. Touboul J, Kharif C. Nonlinear evolution of the modulational instability under weak forcing and damping. *Natural Hazards and Earth System Science* 2010; **10**:2589–2597. DOI: 10.5194/nhess-10-2589-2010.
7. Ducrozet G, Bonnefoy F, Le Touzé D, Ferrant P. 3-D HOS simulations of extreme waves in open seas. *Natural Hazards and Earth System Sciences* 2007; **7**(1):109–122.

8. Bingham H, Zhang H. On the accuracy of finite-difference solutions for nonlinear water waves. *Journal of Engineering Mathematics* 2007; **58**:211–228.

9. Engsig-Karup A, Bingham H, Lindberg O. An efficient flexible-order model for 3D nonlinear water waves. *Journal of Computational Physics* 2009; **228**:2100–2118.

10. Asaithambi NS. Computation of free-surface flows. *Journal of Computational Physics* 1987; **73**:380–394.

11. DeSilva S, Guenther R, Hudspeth R. Irregular points in 2-D free surface flows with surface tension for the wavemaker boundary value problem. *Applied Ocean Research* 1996; **18**:293–302.

12. Li B, Fleming C. A three dimensional multigrid model for fully nonlinear water waves. *Coastal Engineering* 1997; **30**:235–258.

13. Zakharov V. Stability of periodic waves of finite amplitude on the surface of a deep fluid. *Journal of Applied Mechanics and Technical Physics* 1968; **9**:190–194.

14. Engsig-Karup A. Efficient low-storage solution of unsteady fully nonlinear water waves using a defect correction method. *Journal of Scientific Computing* 2011. Submitted.

15. Ducrozet G, Bonnefoy F, Le Touzé D, Ferrant P. Implementation and validation of nonlinear wave maker models in a HOS numerical wave tank. *International Journal of Offshore and Polar Engineering* 2007; **16**:161–167.

16. Ducrozet G, Bonnefoy F, Le Touzé D, Ferrant P. A modified high-order spectral method for wavemaker modeling in a numerical wave tank. *European Journal of Mechanics - B/Fluids* 2011. Accepted for publication.

17. Bonnefoy F, Ducrozet G, Le Touzé D, Ferrant P. Time-domain simulation of nonlinear water waves using spectral methods. In *Advances in Numerical Simulation of Nonlinear Water Waves*, Vol. 11, Advances in Coastal and Ocean Engineering. World Scientific, 2009; 129–164.

18. Rienecker M, Fenton J. A Fourier approximation method for steady water waves. *Journal of Fluid Mechanics* 1981; **104**:119–137.

19. Fornberg B. *A practical guide to pseudospectral methods*. Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, 1995.

20. Fructus D, Clamond D, Grue J, Kristiansen Ø. An efficient model for three-dimensional surface wave simulations. Part I: Free space problems. *Journal of Computational Physics* 2005; **205**:665–685.