# A non-linear wave decomposition model for efficient wave-structure interaction.

## Part A: Formulation, validations and analysis

Guillaume Ducrozeta, Allan P. Engsig-Karupb, Harry B. Binghamc, Pierre Ferranta

<sup>a</sup>LHEEA Lab., École Centrale de Nantes, LUNAM Université, Nantes, France

<sup>b</sup>Dept. of Informatics and Mathematical Modelling, Technical University of Denmark, Kgs. Lyngby, Denmark

<sup>c</sup>Dept. of Mechanical Engineering, Technical University of Denmark, Kgs. Lyngby, Denmark

### Abstract

This paper deals with the development of an enhanced model for solving wave-wave and wave-structure interaction problems. We describe the application of a nonlinear splitting method originally suggested by Di Mascio et al. [1], to the high-order finite difference model developed by Bingham et al. [2] and extended by Engsig-Karup et al. [3, 4]. The enhanced strategy is based on splitting all solution variables into incident and scattered fields, where the incident field is assumed to be known and only the scattered field needs to be computed by the numerical model. Although this splitting technique has been applied to both potential flow and Navier-Stokes solvers in the past, it has not been thoroughly described and analysed, nor has it been presented in widely read journals. Here we describe the method in detail and carefully analyse its performance using several 2D linear and non-linear test cases. In particular, we consider the extreme case of non-linear waves up to the point of breaking reflecting from a vertical wall; and conclude that no limitations
are imposed by adopting this splitting. The advantages of this strategy in
terms of robustness, accuracy and efficiency are also demonstrated by com-
parison with the more common strategy of solving the incident and scattered
fields together.

Keywords: High-Order Finite Differences, OceanWave3D, Wave-structure interaction, Non-Linear Decomposition, Standing waves, Offshore Engineering

# 1. Introduction

The unsteady problem of the interaction between waves and marine struc-
tures is a subject of major interest considering the economic stakes associ-
ated with it. The design of these structures (e.g. ships, offshore platforms,
breakwaters) in coastal and open ocean regions requires very specific tools
to evaluate the expected hydrodynamic forces. Experimental tests on model
scale are considered the most reliable way to obtain such data. However,
due to developments made within the scientific community in recent years,
numerical models are now able to achieve similar levels of reliability.

To obtain the most accurate numerical predictions, the properties of non-
linear phenomena involved may be taken into account. These include: wave
dispersion, wave-wave interaction, diffraction/radiation and interaction with
the bathymetry. When dealing with the development of a numerical method,
one of the key points to take into account for practical use is the computa-
tional work effort. Indeed, for many existing methods this may prohibit the
desired applications, for instance, of shape optimization or parametric stud-
ies. To support cost-efficient analysis with sufficient accuracy there is an
important need for efficient computational models.

Several numerical solutions for coastal and marine wave-structure inter-
action problems have been developed. These rely primarily on potential flow
theory with formulations ranging from linear frequency domain (with pos-
sible developments to second and third-order in non-linearity, e.g. [5]) up
to fully non-linear time-domain [6]. Reynolds-Averaged Navier-Stokes Equa-
tions (RANSE) solvers have also appeared (e.g. [7]). However, the compu-
tational effort required for RANSE methods still imposes severe limitations
on domain size and/or resolution. Thus, we focus here on fully non-linear
time-domain potential flow solvers. The pioneering work of Longuet-Higgins
and Cokelet [8] led to Boundary Integral Equations Methods which allow
a projection of the 3D problem onto the boundary surface of the fluid do-
main. Boundary Element Methods (BEM) have been widely used to solve
this fully non-linear problem in the time-domain (see e.g. recent applications
in [9, 10]). However, efficient BEM solvers are yet to be demonstrated, even
though the latest developments dealing with High-Order BEMs coupled with
Fast Multipole Algorithm [11] are promising. An alternative is to solve the
complete 3D Laplace problem. Different solution strategies may be applied
including: finite element methods [12, 13] and finite difference methods [14].

Recent developments (Bingham & Zhang [2], Engsig-Karup *et al.* [3, 4]
and Glimberg *et al.* [15]) demonstrate the attraction of the direct method
of solving the Laplace problem using high-order finite differences in terms
of both accuracy and efficiency. This method, named OceanWave3D, cap-
tures all the features of wave-wave, wave-bottom and wave-structure inter-
action. However, when studying the interaction between waves and marine
structures, for optimum efficiency it is attractive to minimize the size of the
computational domain around the structure of interest. To this end, an in-
teresting approach was proposed by Di Mascio et al. [1]. They proposed a
non-linear decomposition of the wave-structure interaction solution into an
incident and a scattered wave field. The incident field represents the undis-
turbed wave, which is known explicitly, either from theory or specialized,
fast, and high-accuracy computations; while the scattered field takes into
account the presence of the structure and is solved by the general purpose
numerical model. This approach has been extensively validated for potential
flow problems by Ferrant [16, 17] and also applied to viscous RANSE models
by Ferrant et al. [18] and Luquet et al. [19, 20]. These references are how-
ever not widely read outside of a small Naval Hydrodynamics community,
and they lack a detailed investigation of the features and possible limitations
of the approach. One goal of the present paper is to provide such a detailed
investigation, while at the same time exposing the idea to a wider community
of researchers interested in wave-structure interaction.

In the remainder of this paper, we begin with a brief review of the general
potential formulation as well as the original OceanWave3D model. The non-
linear decomposition is then detailed, where in our case the incident part
is explicitly given by a dedicated wave model, while the scattered part is
solved by the finite difference method OceanWave3D. A challenging feature
of this splitting is that the incident wave field needs to be evaluated up to
the total free-surface height, *i.e.* at points which are above (and sometimes
very far above) the incident wave elevation itself. Fortunately, potential flow
solutions can generally be continued analytically into this region; but it raises
the question of what impact this feature has on the robustness and stability
of the procedure compared to solving for both fields together.

The theory is followed by validation results in 2D. Two configurations are
studied: nonlinear, free wave propagation; and reflection of nonlinear waves
from a vertical wall. The latter case is designed to test the above mentioned
question of robustness and stability by pushing the resulting standing wave
up to the point of breaking, and comparing the results to both experimental
measurements and to numerical calculations where both fields are solved
together. From these test cases, we conclude that the splitting does not
impose any limitations or constraints on the numerical model.

Finally, the efficiency of the non-linear decomposition method is studied
and compared to the original OceanWave3D model. The performance bene-
fits as well as possible disadvantages to using this approach are pointed out
using two generic test cases: reflection of waves from a vertical wall and wave
shoaling on a beach.

## 2. Numerical method

In this section we present the model developed for enhanced treatment
of wave-structure interactions. Firstly, the high-order finite difference model
OceanWave3D [2, 3] is briefly detailed including a short review of the fully
non-linear potential flow solution process. Then, the non-linear decompo-
sition approach is detailed, highlighting its main features and specificities.
The gentle smoothing procedure used to prevent instability when spurious
sawtooth instabilities appear during simulations of steep events is finally dis-
cussed.

## 2.1. General framework - The Ocean Wave3D model

We choose a Cartesian coordinate system with an origin $O$. 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 assume that the fluid simulated is incompressible and inviscid and that the flow is irrotational (i.e. potential flow formalism). With these assumptions, the velocity $\mathbf{V}$ derives from a velocity potential $\phi$. The continuity equation thus becomes the Laplace equation in the fluid domain $D$: $\Delta\phi = 0$.

The free surface elevation is assumed to be a single-valued and thus described by $z = \zeta(\mathbf{x}, t)$. At the same time we introduce, following Zakharov [21], the free surface velocity potential $\tilde{\phi}(\mathbf{x}, t) = \phi(\mathbf{x}, \zeta, t)$. This way, the fully non-linear time-dependent free surface boundary conditions can be written in terms of surface quantities

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

$$ \zeta_t = (1 + |\nabla \zeta|^2) \phi_z - \nabla \tilde{\phi} \cdot \nabla \zeta \quad (2) $$

on $z = \zeta(\mathbf{x}, t)$; where partial differentiation is indicated when the independent variables appear as sub-scripts and $\nabla$ stands for the horizontal gradient. It appears in these two conditions that one only need to evaluate the vertical velocity at the free surface $\tilde{w}(\mathbf{x}, t) = \phi_z(\mathbf{x}, \zeta, t)$ in order to advance in time the quantities of interest $(\zeta, \tilde{\phi})$.

The OceanWave3D model solves the Laplace problem (i.e. Laplace equation throughout the depth of the fluid along with boundary conditions - including bottom whose profile can be chosen arbitrarily) directly to obtain
this vertical velocity. A non-conformal mapping ($\sigma$-transform) allows the transformation of the physical domain (which evolves during time) onto a fixed computational domain. The obtained Laplace problem is then solved by means of finite-difference method. Arbitrary high-order schemes are used with possibly stretched grid in horizontal or vertical directions. These developments are found to be greatly advantageous relative to classical second-order schemes on even grids. The obtained linear system is solved with low storage and high efficiency with a three-recurrence multigrid preconditioned Defect Correction Method [4]. More details can be found in the cited references.

## 2.2. Non-linear decomposition

Efficient numerical computation of wave-body interactions is still a challenge for fully non-linear potential flow methods (such as OceanWave3D). Thus in order to minimize the size of the required computational domain, we adopt the splitting technique originally proposed by Di Mascio et al. [1]. All the unknowns of the problem (velocity potential and free-surface elevation) are split into two parts which are referred as the incident and the scattered components. The incident terms describe the undisturbed incoming wave and the scattered terms take care of the influence of the structure. The incident component is assumed to be known explicitly while the scattered component is the new unknown to be solved for by the OceanWave3D model. The incident waves are described explicitly using a non-linear potential flow model and will act as forcing terms. In the following results, non-linear regular wave trains are studied and given by an algorithm based on the stream function theory of Fenton [22] (or Airy waves when solving in
a linear framework). This approach has been extensively validated by Ferrant [16, 17] in the XWAVE code which used this splitting technique with a MEL-BEM method. A successful comparison between this XWAVE code and the modified OceanWave3D code has been presented in Ducrozet et al. [23]. Furthermore, the approach has been successfully extended to viscous flow calculations with the so-called SWENSE method (Ferrant et al. [18], Luquet et al. [19, 20]).

We apply the following change of variable to the total free surface elevation $ζ$ and the velocity potential $φ$:

$$ \zeta = \zeta^I + \zeta^S \qquad (3a) $$

$$ \phi = \phi^I + \phi^S \qquad -h \le z \le \zeta \qquad (3b) $$

The incident waves (described by $ζ^I$ and $φ^I$) are known potential flow solutions (thus $Δφ^I = 0$ in the fluid domain $D$ where $Δ$ is the 3D Laplacian). The modified time-dependent free surface boundary conditions now read:

$$ \zeta_t^S = -\zeta_t^I - \nabla (\zeta^I + \zeta^S) \cdot \nabla (\tilde{\phi}^I + \tilde{\phi}^S) \\ + (\tilde{w}^I + \tilde{w}^S) (1 + \nabla (\zeta^I + \zeta^S) \cdot \nabla (\zeta^I + \zeta^S)) \qquad (4a) $$

$$ \tilde{\phi}_t^S = -\tilde{\phi}_t^I - g(\zeta^I + \zeta^S) - \frac{1}{2} \nabla (\tilde{\phi}^I + \tilde{\phi}^S) \cdot \nabla (\tilde{\phi}^I + \tilde{\phi}^S) \\ + \frac{1}{2} (\tilde{w}^I + \tilde{w}^S)^2 (1 + \nabla (\zeta^I + \zeta^S) \cdot \nabla (\zeta^I + \zeta^S)) \qquad (4b) $$

expressed on $z = \zeta^I + \zeta^S$, the total free surface elevation.

The Laplace problem becomes

$$ \nabla^2 \phi^S + \phi_{zz}^S = 0, \qquad -h \le z < \zeta^I + \zeta^S \qquad (5a) $$

$$ \phi^S = \tilde{\phi}^S, \qquad z = \zeta^I + \zeta^S \qquad (5b) $$

$$ \vec{n} \cdot \vec{\nabla} (\phi^I + \phi^S) = 0, \qquad (\mathbf{x}, z) \in \partial\Omega \qquad (5c) $$
where $\partial\Omega$ represents the solid boundaries of the domain $\vec{\mathbf{n}}$ is the unit normal vector to the corresponding boundary and $\vec{\nabla}$ the 3D gradient operator. Notice that this makes the kinematic boundary conditions on $\partial\Omega$ inhomogeneous in general. Otherwise, the Laplace problem is almost unchanged compared to the original version. However, it is solved using the $\sigma$-transformation which will introduce several changes in the solution. Inserting the decomposition into the transformation gives

$$ \sigma(\mathbf{x}, z, t) = \frac{z + h(\mathbf{x})}{\zeta^I(\mathbf{x}, t) + \zeta^S(\mathbf{x}, t) + h(\mathbf{x})} = \frac{z + h(\mathbf{x})}{d(\mathbf{x}, t)} $$

with $d(\mathbf{x}, t) = \zeta^I(\mathbf{x}, t) + \zeta^S(\mathbf{x}, t) + h(\mathbf{x})$. The Laplace problem is thus transformed to (introducing for clarity $\Phi(\mathbf{x}, \sigma, t) = \phi(\mathbf{x}, z, t)$)

$$ \Phi^S = \tilde{\phi}^S, \quad \sigma = 1 \quad (6a) $$

$$ \nabla^2 \Phi^S + \nabla^2 \sigma \Phi^S_{\sigma} + 2 \nabla \sigma \cdot \nabla \Phi^S_{\sigma} + (\nabla \sigma \cdot \nabla \sigma + \sigma_z^2) \Phi^S_{\sigma\sigma} = 0, \quad 0 \le \sigma < 1 \quad (6b) $$

$$ (\sigma_z + \nabla h \cdot \nabla \sigma) \cdot (\Phi^I_{\sigma} + \Phi^S_{\sigma}) + \nabla h \cdot \nabla (\Phi^I + \Phi^S) = 0, \quad \sigma = 0 \quad (6c) $$

The derivatives of $\sigma$ are expressed with the change of variable as

$$ \nabla \sigma = (1 - \sigma) \frac{\nabla h}{d} - \sigma \frac{\nabla (\zeta^I + \zeta^S)}{d} \quad (7) $$

$$ \nabla^2 \sigma = \frac{1-\sigma}{d} \left( \nabla^2 h - \frac{\nabla h \cdot \nabla h}{d} \right) - \frac{\sigma}{d} \left( \nabla^2 (\zeta^I + \zeta^S) - \frac{\nabla (\zeta^I + \zeta^S) \cdot \nabla (\zeta^I + \zeta^S)}{d} \right) - \frac{1-2\sigma}{d^2} \nabla h \cdot \nabla (\zeta^I + \zeta^S) - \frac{\nabla \sigma}{d} \cdot (\nabla h + \nabla (\zeta^I + \zeta^S)) \quad (8) $$

$$ \sigma_z = \frac{1}{d} \quad (9) $$
The discretization of the model has been studied in detail. The Laplace problem (Eqs. 6) as well as free surface boundary conditions (Eqs. 4) make use of incident quantities and their spatial derivatives. Then, the discrete model for solving the problem of interest can be developed in two ways: i) use the analytical derivatives directly as source contributions or ii) sum the computed solution for the scattered field with the incident field and then do all approximations discretely. The approximate solution is more straightforward in terms of implementation and requires less computational effort to generate the incident wave field. However, this extra computational effort is expected to be negligible in 3D computations where the resolution of scattered field will be predominant. At the same time, this approach leads to accuracy losses. Indeed, different tests conducted showed that both approaches converge to a similar behavior (at large number of points), but the analytical derivatives are far more accurate for a small number of points. Consequently, we choose to use the analytical description of all incident quantities in our implementation of the non-linear decomposition.

## 2.3. The main features of the decomposition

The proposed non-linear decomposition approach presents several noteworthy features. They are described hereafter in detail to point out clearly the advantages and possible limitations that are expected from such an approach. We also refer to the different validations proposed in Sec. 3 for practical applications of some of these features.

First of all, the incident wave field is given explicitly by a dedicated wave model which is assumed to be accurate and very fast (compared to the overall
solution). As noted above, a stream function model is used in the following
to generate non-linear regular waves however, the choice of incident wave
model is arbitrary. Consequently, the only limitation on the incident wave
field is the capability of the chosen wave model. The SWENSE approach [20]
has for instance been coupled to a High-Order Spectral (HOS) method [24]
to define complex sea states (including irregular or focused 2D and 3D wave
fields).

Concerning the scattered component which is actually solved by the mod-
ified OceanWave3D solver (i.e. scattered potential φS and free surface eleva-
tion ζS), it seems important to detail what is included in this part. Thanks to
the non-linear decomposition, the incident part is given explicitly in the set
of equations solved for the scattered part. Consequently, physical processes
which are not included in the incident part are solved for by the scattered
component. Therefore, this scattered part is not only the so-called scattered
(or diffracted) waves by the object of interest but also all interactions be-
tween these waves and the incident components. Thus, it is advantageous to
choose an incident wave field which solves as much of the problem as pos-
sible. In particular, it is worthwhile to use a fully non-linear model for its
description.

Indeed, if a linear incident wave field is chosen instead, all non-linearities
which are not included by incident wave model will be solved by the "scat-
tered" part. Similarly, if the incident wave field is designed at a given water
depth and one wants to model the refraction of waves by an uneven bottom,
this refraction process will be actually solved by the scattered part (see Sec.
4.2). On the other hand, if all physical processes are included by the incident part, the corresponding scattered part is zero and requires no CPU cost (see Sec. 3.1). Consequently, the more of the physics that are resolved by the incident part, the more efficient the methodology becomes.

The main advantages which are expected may be listed from the previous features. Firstly, the computational effort of the incident wave field will be negligible (especially in 3D) compared to the effective solution of the wave-structure interactions. Furthermore, we note that with an appropriate choice of incident wave model, any kind of wave can be defined and generated. Also, there is no need to propagate the incident waves from the computational domain boundary to the structure with this approach. Indeed, the incident waves exist everywhere in the domain at all times. We note however, that the original OceanWave3D model offers already this possibility via a special use of relaxation zones.

In addition, the incident waves are not reflected from the far-field boundaries. Thanks to the non-linear decomposition, it is then also possible to reduce the problem size for 3D scattering problems. Since the incident flow is known everywhere with no reflection from exterior boundaries, the grid may be highly refined around the structure, and dramatically stretched away from it, leading to more efficient use of the computational unknowns. An accurate solution of the scattered field is obtained close to the structure, but the domain away from the structure only functions to absorb the scattered waves. This results in a reduced problem size and a corresponding reduction of the computational effort.
Consider the free surface boundary conditions (4). From these expres-
sions we can see that the incident field and its derivatives ( φ<sup>I</sup>, φ<sub>t</sub><sup>I</sup>, w<sup>I</sup>, ... )
have to be evaluated at the position of the total wave elevation ζ = ζ<sup>I</sup> + ζ<sup>S</sup>.
Thus, these incident quantities have to be evaluated at points above the
incident free surface elevation (as sketched on Fig. 1). Consequently, the
dynamics of the incident wave has to be computed outside its domain of def-
inition. The accuracy of this evaluation is in fact not critical since everything
which is not solved by the incident part will be solved by the scattered part.
Theoretically it should be possible to specify an incident wave field only in
the incident fluid domain and set it to zero outside, for example. However,
numerically it appears to be necessary to have a continuous solution for the
incident potential above the incident free surface elevation to ensure that dis-
continuities are not introduced in the scattered potential. This also leads to
a more efficient and stable solution. Thus the ideal incident wave field, and
its derivatives, should be continuouly defined through, and above the free-
surface. Fortunately, potential flow based wave theories (such as the stream
function theory used here or High-Order Spectral method [24]) generally al-
low for an analytic continuation of the incident potential above the reference
wave elevation. Thus, this continuity requirement is not a constraint in prac-
tice. However, it is natural to question whether errors introduced as a result
of this extrapolation (which will be corrected by the scattered part) may be
a source of instabilities or inaccuracies. This has motivated the choice of
the validation test case in Sec. 3.2 with the reflection of waves by a vertical
wall; which is assumed to be the most extreme situation we might expect to
encounter in this respect.
<figure><img src="image_1.png" /><figcaption>Figure 1: Illustration of a possible situation showing how the incident wave may need to be evaluated at points far above its own free-surface definition</figcaption></figure>

Finally, the initialisation of such an approach has to be treated carefully. In order to ensure a stable solution, a time ramp is imposed at the beginning of the computation on the non-homogeneous Neumann boundary conditions (5c). These are indeed imposed on the combined incident and scattered field and require special care during the start-up. Thus, the problem solved is going from a wave propagating freely inside the domain at the beginning of the time ramp up to a wave interacting with the structure at the end: the structure studied slowly *appears* during the time-ramp. This solution strategy has been found to be robust, and effective at preventing numerical instabilities. In the following, the length of this time ramp is usually fixed to 5 *T* with *T* the characteristic wave period with a solution independent to time ramp after 10 *T*.
## 2.4. The smoothing procedure

Here we detail the smoothing procedure used during the nonlinear computations which involve very steep waves (close to wave breaking) and possible spurious sawtooth instabilities. Inside the domain (i.e. away from boundaries), we use 1D Savitzky Golay (SG) filters [25] applied along each direction successively for 2D surface smoothing (i.e. in 3D cases). Kuo et al. [26] demonstrated the equivalence of sequential 1D smoothing and 2D smoothing techniques. This is briefly reviewed next. Then we study the efficient treatment of boundaries during the smoothing. The classical off-centered SG-filters are analysed and switched to optimized filters from Berland et al. [27] to obtain better stability properties.

Note that the proposed smoothing is applied only in extreme conditions and occasionally to ensure that it does not influence significantly the total energy. Typically, it is applied every wave period to every tenth of a wave period, depending on the resolution and the test case.

### 2.4.1. Centered Savitzky-Golay filter

In this paragraph we review the behavior of 1D centered Savitzky-Golay filters [25]. This filter will be used, if needed, inside the domain away from boundaries. Figure 2 presents the behavior of the amplification factor with respect to the number of points per wavelength describing a wave. The amplification factor $A$ is the ratio between the original wave amplitude and the filtered one. A value of 1 means that the wave is not modified by the filter while a value 0 means a complete absorption of the corresponding wave. In this figure, bottom graphs present the value $\log(1-A)$ to have a more accurate view on the strength of the filter.
<figure><img src="image_2.png" /><figcaption>Figure 2: Amplification portraits for SG(13,*p*) filters (left) and SG(*N*,4) filters (right)</figcaption></figure>

We will refer to these filters as SG(*N*,*p*) filters where *N* is the number of grid points used and *p* is the polynomial order of the interpolation. The top left part of Fig. 2 shows the dependence of the amplification factor for SG(13,*p*) filters using different orders of polynomial *p* (from 2<sup>nd</sup> to 10<sup>th</sup> order). The 10<sup>th</sup> order filter is very smooth and only slightly damps the shortest waves with resolution below 5-6 points per wavelength while the 2<sup>nd</sup> order one is very strong and applies significant damping to waves resolved with around 30 points per wavelength. The top right part of Fig. 2 gives the
amplification factors for SG(N,4) filters with N from 7 to 15 points. More
points in the filter leads to stronger filtering. Thus, the two parameters
can be modified to obtain a desired behavior of the filter. Typically, we
employ the mild filter SG(13,10) which for most purposes is sufficient for
stabilization. It is applied at a given frequency depending on the intensity of
filtering needed. The bottom part of the figure with log diagram assess that
this filter is very light, even for long time simulations where it can be applied
numerous times. In all cases it is important to monitor the conservation of
energy to ensure that the physics is not contaminated by the filtering.

2.4.2. Boundary treatment

Towards the boundaries of the domain, off-centered schemes are required.
In this Section, we analyse the properties of off-centered Savitsky-Golay fil-
ters and show that none of them are suitable very close to the boundary
where they act to amplify rather than filter. Thus, we adopt instead the
selective filters of Berland *et al.* [27].

We study the behavior of the off-centered SG filters and particularly the
one which is completely one-sided (i.e. the 1st point of the domain - i.e.
boundary point - with neighbours all on one side). It is then possible to draw
amplification portraits similarly to previous figures. To draw such a portrait,
one applies the same off-centered scheme everywhere in the domain. This is
actually not exactly equivalent to the current filter which has a progressive
off-centering once one reaches the boundary and a centered scheme far from
the boundaries. However, it allows us to characterize its behavior at the
boundary position. This is done in Fig. 3.
<figure><img src="image_3.png" /><figcaption>Figure 3: Amplification portraits on 1st point for off-centered filters SG(13,p) (left) and SG(N,4) (right)</figcaption></figure>

Different configurations are depicted in Fig. 3 (different choices of polynomial order and number of points in the filter) and a major issue appears, noticing that this SG off-centered filter actually amplifies some wave components instead of damping them. In addition, high-order polynomial (e.g. SG(13,10)) exhibits no filtering at all on the point of interest. This has driven the choice to use an alternative smoothing procedure.

To overcome these problems, some optimized filters have been developed by Berland et al. [27], called Selective Filters (SF). These are some explicit filters formulated to treat the boundary problems with some specifications and particularly the one of interest, that no amplification occurs at any wavelength. The filtered data $\tilde{f}(x_i)$ is expressed with respect to the unfiltered data $f(x_i)$ as

$$ \tilde{f}(x_i) = f(x_i) - \alpha \sum_{j=-a}^{b} d_j f(x_i + j \Delta x) \quad (10) $$
For now, we set the parameter $\alpha = 1$ to analyze the behavior of the filter. However, it will be possible to adjust it to have a lighter smoothing if needed ($\alpha$ has to be chosen between 0 (no filtering) and 1), but note that the features described hereafter are preserved.

Figure 4 presents the amplification portraits of these filters which we refer to as $SF(a, b)$, with $a$ and $b$ the number of points on the left and on the right respectively. Thus, $SF(0, 3)$ stands for a one sided scheme that can be applied on the left boundary itself.

<figure><img src="image_4.png" /><figcaption>Figure 4: Amplification portraits of Selective Filters from Berland et al. [27]: normal (left) and log scale (right)</figcaption></figure>

This figure (left) indicates that these filters do not exhibit any amplification. They also appear to be really light, with an amount of dissipation of the order of the $SG(13, 10)$ filter. It is then possible to choose a set of filters at every off-centered position to have always almost the same amount of dissipation everywhere in the domain (when using smoothing).
In order to assess more explicitly the difference between the filters, Fig. 4 (right) plots the amplification portraits on a log-log scale. It is more clear here that the correct behaviour is achieved at the expense of some accuracy with long waves. However, this is not expected to be a problem, since these filters are only applied to a very small number of grid points near the boundaries.

Several tests have been performed to deduce correct parametrization of these boundary-filters. In the following, we typically use a given $\alpha$ value which is taken lighter on the first two points from boundary ($\alpha_1 = \alpha_2 = \alpha/10$) when smoothing is applied (same approach used in Berland et al. [27]). $\alpha_i$ stands for the value of the parameter on the $i^{th}$ point ($i = 1$ is the boundary).

### 3. Validation

In this section, the first validation results in 2D domains (i.e. one horizontal direction) are given to demonstrate the feasibility of decomposition presented in the previous section. First, we study the case of a non-linear wave propagating without any obstacles inside the domain. Then, we investigate the case of a mildly non-linear wave reflecting off a wall. To ensure the robustness of the non-linear decomposition, an extreme case is simulated in this reflection configuration. The limit of stability of the standing wave appearing close to the wall is discussed and compared with experimental results as well as a dedicated numerical model. All these validations take the incident wave to a non-linear regular wave from stream-function theory [22].

#### 3.1. *Non-linear propagation*

We consider the simple case of a non-linear wave freely propagating without reflection from the side walls of the domain. The decomposition approach
is used with the solution of the problem being specified as a mildly non-linear
incident wave with $H/\lambda = 0.05$ and $kh = 2\pi$, where $\lambda$ is the wave length $H$
the wave height, and $k = 2\pi/\lambda$ the wave number. The results are shown in
Fig. 6 computed using $N_x = 16$ points per wavelength and $N_z = 10$ points
in the vertical. A stretched grid is used in the vertical direction as sketched
in Fig. 5. The iterative DC Method is used with a tolerance level of $10^{-13}$
and a wavelength $\lambda = 1$ is chosen.

<figure><img src="image_5.png" /><figcaption>Figure 5: Example of mesh used with OceanWave3D</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6: Propagation of mildly non-linear wave $H/\lambda = 0.05$. Free surface elevation: incident $\zeta^I$ and scattered $\zeta^S$ after 25 periods of propagation</figcaption></figure>

In this figure, the left part shows the incident free surface elevation $\zeta^I$ after
25 periods of propagation. The right part gives the corresponding scattered
free surface elevation $\zeta^S$ which is actually the part solved directly by the
enhanced OceanWave3D model. The solution of the problem of interest is the sum of incident and scattered parts: $\zeta = \zeta^I + \zeta^S$.

Clearly the scattered wave field is close to zero in the whole domain ($<$ $10^{-11}$). The amplitude of this diffracted field will be strongly correlated with the discretization parameters of the numerical method (e.g. tolerance for iterative solvers as well as truncation errors resulting from discretization of governing equations). Indeed, in the case of a freely propagating wave field, we know (this computation enables a validation) that the incident wave field strictly satisfies the non-linear system Eqn. (6) to be solved by the direct method and is consequently a proper solution to the problem. Therefore, the coupled problem solved by direct method will just check that the incident wave field satisfies this system to the level of truncation errors. Thus, the scattered field created corresponds only to some minor errors in this verification (correlated with the tolerance of the DC method). This appears also clearly when one investigates the number of iterations used to solve the linear system. The latter is always equal to 1, independently of the space and time discretization used or the tolerance chosen, indicating that incident wave field is a solution to the problem considered. We return to a discussion of computational effort in Sec. 4.1.3.

In contrast, if one uses in such a non-linear simulation an incorrect incident wave field (i.e. a sea state that does not satisfy the system of equations (6)) defined as a linear one for instance, then non-linearities will appear in the scattered part. Similarly, with the incident field defined here at a fixed water depth (from stream function theory), if we study the case of a variable
bathymetry the influence of the bottom will be described in the scattered
wave field as discussed later (Sec. 4.2). Another point of view is that if inci-
dent field is defined to be zero, then the solver will correspond to the original
solver since all effects will be resolved in the scattered variables. Thus, this
extension simply adds to the functionality of the original solver.

3.2. *Non-linear reflection on vertical wall*

The reflection of a wave on a vertical wall is an interesting case, partic-
ularly to test both the robustness and the accuracy of the decomposition
approach. Indeed, as seen in Sec. 2.3, we have to evaluate the incident
quantities at the position of the total free surface elevation. This introduces
some evaluations outside its domain of definition and is a possible source of
errors or instability during the solution. In this respect, the reflection of a
wave from a wall is an extreme situation where the incident field must be
evaluated at a position which is at least twice the incident wave amplitude.
Furthermore, strong non-linear interactions are expected between the inci-
dent and reflected waves. These interactions must also accounted for in the
scattered part. They are expected to be described with complex frequency
content and consequently the scattered part will be more complicated than
just a simple propagating wave in the opposite direction (w.r.t. the incident
one). Thus, it seems reasonable to think that if the method is capable of
handling this particular case with high incident wave amplitudes, it will be
able to treat any non-breaking wave-structure interaction problems that we
will have to deal with.

Figure 7 shows the results of a calculation using the same mildly non-
linear wave as in the previous test case: $H/\lambda = 0.05$ and $kh = 2\pi$ (i.e.
deep water). The computational domain is of length $L_x = 6\lambda$ with the waves propagating from right to left and reflecting from the left boundary. An absorbing zone of length $4\lambda$ is located on the right part of the domain to absorb the scattered waves and prevent them from reflecting on the right wall. The spatial discretizations are $N_x = 32$ points per wavelength and $N_z = 10$ points in the horizontal and vertical direction respectively. A wavelength $\lambda = 1m$ is chosen for this figure.

Results are presented similarly to Fig. 6. Each row presenting two columns gives, at a given time step (indicated on the left side), the incident wave $\zeta^I$ and the corresponding scattered component $\zeta^S$ on the left and right part of the figure respectively. Steady-state results are presented for 4 time steps selected during one incident wave period $T$. The bottom row presents the total free surface elevation $\zeta = \zeta^I + \zeta^S$ with several time steps superimposed (corresponding to the different colors) chosen also during one wave period. This representation is intended to point out the creation of a steep non-linear standing wave close to the left wall. Reminding that the absorption is applied to the scattered part in the range $x \in [2; 6]$, the standing wave is actually limited to the region close to the left wall ($x \in [0; 2]$). In this figure we observe, as expected, that the scattered wave field is complex but also of high amplitude (larger than the incident one). Non-linearities appear to be essential noticing that in the linear case the scattered wave will be the same than the incident one, with opposite direction. In contrast, this non-linear test case induces a scattered part $\zeta^S$ very different from the incident wave. However, the resulting wave field $\zeta$ solution of the problem of an incident regular wave reflecting on a vertical wall is, as expected, a
<figure><img src="image_7.png" /><figcaption>Figure 7: Non-linear reflection on a vertical wall $H/λ = 0.05$. Free surface elevation: incident $ζ^I$, scattered $ζ^S$ and total for different time-steps. Bottom figure captures free surface elevation at different time steps.</figcaption></figure>
standing wave. We can observe on Fig. 7 that non-linear effects result in a standing wave whose wavelength is reduced compared to the incident wave: $\lambda_{stand.} = 0.95\lambda$. The wave height is also amplified compared to the linear case: $H_{stand.}/\lambda = 0.107$, i.e. an effective steepness $H_{stand.}/\lambda_{stand.} = 0.113$, i.e. 13% larger than the linear solution.

The difference of wavelength between the incident wave $\zeta^I$ and the effective solution of the wave-structure interaction $\zeta$ imposes a complex scattered wave field with spatially varying amplitude as observed in Fig. 7. Furthermore, if we take a look at the phasing of the two components arising from non-linear decomposition ($\zeta^I$ and $\zeta^S$) at the position of the wall ($x=0$), they are clearly out of phase. In addition, during this test case, the maximum of the scattered component $\zeta^S$ is observed when the incident wave is almost equal to zero on the wall ($t = t_0 + T/2$) in Fig. 7. This phase shift explains the high-amplitude consequently observed in the scattered wave field. The phasing of the two components is studied specifically in Fig. 8 which presents the temporal evolution of the free surface elevation on the vertical wall where the reflection occurs.

A time-ramp is applied at the beginning of the simulation during $10T$ on the non-homogeneous Neumann boundary conditions. This ensures a stable solution for all cases treated with this configuration. The figure denotes clearly that when the steady-state is reached (after $t \approx 20T$), the incident and the scattered wave fields are phase-shifted (according to what we observed on free surface elevation plots Fig. 7 and opposite to the linear case). This induces some high amplitudes of the scattered field (compared to the incident one) to obtain the solution of the problem which is consequently also
<figure><img src="image_8.png" /><figcaption>Figure 8: Free surface elevation on the vertical left wall during time until steady state is reached: incident $\zeta^I$, scattered $\zeta^S$ and total $\zeta$</figcaption></figure>
phase shifted compared to the incident wave field. This test-case thus turns
out to be even more extreme than originally envisioned (with respect to the
extrapolation of incident quantities). Indeed, these have to be evaluated far
above the incident wave elevation (for instance $t = 27.2$ s, $\zeta^I = 0.01$ and
$\zeta^I + \zeta^S = 0.065$). However, this phenomenon does not induce any stability
or robustness issues compared to the behavior of the original OceanWave3D
model. Various tests have been made to ensure this conclusion, including
increasing the discretizations in the two space-directions as well as the local
refinement near the wall where the reflection occurs.

Then, it seems interesting to study the influence of non-linearities on the
phase shift observed previously. This is done in Fig. 9 which presents the
phase shift between the two components at the location of the wall. This is
studied as a function of the steepness of the incident wave field $H/\lambda$.

<figure><img src="image_9.png" /><figcaption>Figure 9: Evolution of the phase shift at the wall as a function of incident wave steepness $H/\lambda$</figcaption></figure>
We observe, as expected, that for low incident steepness the phase shift tends to zero as predicted by linear theory. The phase-shift then increases with increasing incident wave non-linearity. It may be seen as the result of the non-linear interaction between incident and reflected waves. This is similar to what is observed during the interaction of a solitary wave with a vertical wall. It has been shown that a phase shift, dependant on spatial position, is observed between the incident and reflected waves. This is typically what is seen here in the case of deep water regular waves. We refer to Su and Mirie [28] or Fenton and Rienecker [29] for details on the solitary-wave configuration.

The results with $H/λ = 0.05$ have been compared to the same case studied with the original OceanWave3D model. The numerical parameters are $N_x = 16$ points per wavelength and $N_z = 10$ points in horizontal and vertical directions respectively with a Courant number $C_r = 0.64$. A relaxation zone on the right side of the domain is used to generate waves with the original model [2]. The two results obtained (with and without the non-linear decomposition) are in excellent agreement. An example of the comparison is given in Fig. 10 at a given time step chosen after a steady state was reached. The original and non-linear decomposition results are superposed on the top part of the figure. The relative difference (in %) is given in the bottom of the figure in the zone of interest ($x \in [0; 2]$) showing that they agree to within 1%. This small discrepancy may be explain by the treatment of reflected waves in the absorption/generation zone (right side of the domain). In this region, the two approaches are not strictly equivalent. Note that the differences in
terms of computational effort are studied in section 4.1.3.

<figure><img src="image_10.png" /><figcaption>Figure 10: Wave profile during reflection on left wall: original vs non-linear decomposition (top) and relative difference in % (bottom)</figcaption></figure>

### 3.3. The stability limit of a standing wave

Having shown the proper solution of a problem with a mildly non-linear incident wave, along with the robustness of the new formulation, it is interesting to determine the limitations of the decomposition in terms of numerical stability. To do this, we now increase the amplitude of the incident wave until the solution becomes unstable. The exact physical stability limit of a standing wave has not been definitively determined yet, but it is approximately $H/\lambda \simeq 0.20$ in deep water cf. [30, 31, 32].

During the different simulations, the choice of numerical parameters is made so that the result is converged and that a steady state is reached.
Furthermore, note that the different results were compared successfully with the original direct method (without decomposition) to ensure the quality of the solution. We find the stability limit to be obtained with an incident wave of steepness $H/\lambda = 0.075$. Figure 11 compares the computed results to experimental measurements by Taylor [33] and to the numerical solution of Okamura [32]. The experimental data correspond to the steepest stable standing wave observed during this experimental campaign. Similarly, the numerical solution of Okamura presents the limiting standing wave whose crest acceleration is 99.98% of gravity. The profiles with a crest in the middle of the domain as well as a trough are given. The numerical simulations are performed here with $N_x = 64$ points per wavelength and $N_z = 10$ points.

<figure><img src="image_11.png" /><figcaption>Figure 11: Standing wave, stability limit: comparison between experiments by Taylor [33] (symbols), numerical solution of Okamura [32] (dashed lines) and numerical simulations with OceanWave3D (straight lines)</figcaption></figure>

Firstly, it may be noted in this figure the good agreement between our
numerical simulation and the experimental results on this extreme case. Troughs appear to be more pronounced in both numerical simulations but we point out the right/left asymmetry of the trough in the experiments, as well as some scatter in the data. Secondly, both numerical results appear to be in good agreement. The model of Okamura [32] is dedicated to the simulation of very steep standing waves in order to deduce their profile and particularly the steepest one which is depicted on the figure. This nice agreement indicates the accuracy of our numerical method. Furthermore, this successful comparison indicates that the proposed decomposition approach is very robust. Indeed, the simulation of this strongly non-linear phenomenon, very close to physical wave breaking, is possible. A light filtering is performed to ensure the stability of the calculation as described in Sec. 2.4. However it is identical to the one required for the original method under the same wave conditions, confirming that the proposed decomposition does not induce any instabilities.

## 4. Efficiency of non-linear decomposition

The introduction of the non-linear decomposition is intended to enhance the efficiency of the original OceanWave3D model. In this section, we described the different performance benefits (computational efficiency, accuracy) one may observe as well as the possible limitations. In the following subsections, we will compare the original OceanWave3D model with the enhanced one. We are ultimately interested in computational effort however we focus on algorithmic efficiency in the immediate comparisons. This algorithmic efficiency will be evaluated thanks to the investigation of the iteration
counts needed in the iterative method (DC or GMRES) for solving the linear
systems within a given tolerance. This choice is motivated by the necessary
resolution requirements to achieve a certain accuracy level overall that may
change between the two approaches. Furthermore, in this 2D configuration,
the additional computational effort to evaluate the incident wave field may
be no longer negligible compared to the overall computation time (which will
happen in 3D). However, notice that it is possible to lower the cost of gener-
ating incident wave fields, e.g. using tabulated data over a wave period for
regular waves. Therefore, we focus mainly on the relative difference in iter-
ation counts between both approaches that will be detailed in the following.

4.1. Performance benefits

In order to assess the performance benefits of the proposed decomposition,
we first have a look to the accuracy of the model to propagate wave fields,
followed by a study of the convergence rate when treating wave-structure
interaction. Then, the efficiency (algorithmic and numerical) is investigated
in those two configurations.

4.1.1. Accuracy

We consider the long-time propagation of non-linear regular waves. In
this test-case, periodic boundary conditions are applied in the horizontal di-
rection. Dispersion errors $\epsilon_\phi$ as well as amplitude ones $\epsilon_a$ are recorded with
a comparison between the original OceanWave3D model and the new model
developed which includes the decomposition. The wave which is considered
is mildly non-linear: $H/\lambda = 0.0955$, i.e. a steepness $ka = 0.30$. The long-
time propagation over a flat bottom (deep water $kh = 2\pi$) of this wave
during 1000 T is computed. Parameters are chosen following the procedure used in [34] which is intended to choose optimal parameters (with respect to CPU timings of the original OceanWave3D model). These parameters (discretizations ($N_x, N_z$), order of the finite differences (FD) scheme $p$, absolute tolerance (Tol.)) as well as corresponding results are reported in Tab. 4.1.1. Comparisons are performed between computed results and analytical results from stream function theory [22]. The average iteration count is obtained using last computed solution as initial guess for each RK-stage.

<table><caption>Table 1: OceanWave3D: Error of phase shift $\epsilon_\phi$ and wave amplitude $\epsilon_a$ after 1000 periods of propagation, $C_r = 0.64$ (i.e. number of time steps going from 12500 ($N_x = 8$) to 400000 ($N_x = 256$)) for a periodic waves ka = 0.3 and kh = $2\pi$</caption><thead><tr><th colspan="4">FD</th><th colspan="3">Original</th><th colspan="3">Enhanced</th><th rowspan="2">Diff.</th></tr><tr><th>$N_x$</th><th>$N_z$</th><th>$p$</th><th>Tol.</th><th>$\epsilon_\phi$</th><th>$\epsilon_a$</th><th>iter.</th><th>$\epsilon_\phi$</th><th>$\epsilon_a$</th><th>iter.</th></tr></thead><tbody><tr><td>8</td><td>7</td><td>6</td><td>$10^{-6}$</td><td>1.3 $10^4$</td><td>7.4 $10^{-1}$</td><td>7.0</td><td>6.8</td><td>1.5 $10^{-2}$</td><td>2.71</td><td>-61%</td></tr><tr><td>16</td><td>7</td><td>6</td><td>$10^{-6}$</td><td>3.7 $10^3$</td><td>2.4 $10^{-1}$</td><td>5.9</td><td>3.5 $10^1$</td><td>4.4 $10^{-2}$</td><td>2.90</td><td>-51%</td></tr><tr><td>32</td><td>9</td><td>6</td><td>$10^{-9}$</td><td>2.3 $10^2$</td><td>2.8 $10^{-2}$</td><td>10.6</td><td>2.7 $10^{-2}$</td><td>3.8 $10^{-5}$</td><td>2.30</td><td>-78%</td></tr><tr><td>32</td><td>9</td><td>8</td><td>$10^{-9}$</td><td>2.3 $10^2$</td><td>2.5 $10^{-2}$</td><td>16.3</td><td>2.1 $10^{-2}$</td><td>4.0 $10^{-5}$</td><td>5.63</td><td>-65%</td></tr><tr><td>64</td><td>17</td><td>8</td><td>$10^{-11}$</td><td>8.1</td><td>2.8 $10^{-3}$</td><td>12.2</td><td>8.8 $10^{-5}$</td><td>2.6 $10^{-6}$</td><td>2.57</td><td>-79%</td></tr><tr><td>128</td><td>33</td><td>8</td><td>$10^{-13}$</td><td>2.6 $10^{-1}$</td><td>5.6 $10^{-5}$</td><td>14.9</td><td>5.6 $10^{-6}$</td><td>2.4 $10^{-6}$</td><td>3.32</td><td>-78%</td></tr><tr><td>256</td><td>65</td><td>8</td><td>$10^{-15}$</td><td>8.6 $10^{-3}$</td><td>3.5 $10^{-6}$</td><td>22.2</td><td>3.7 $10^{-6}$</td><td>2.4 $10^{-6}$</td><td>8.73</td><td>-61%</td></tr></tbody></table>

As expected, both methods converge with the number of points (error on the amplitude of the computed wave as well as its phase shift compared to analytical solution). Furthermore, if one compares both approaches, non-linear decomposition appears more accurate than the original one on this
test case. Particularly, converged results are obtained with smaller number of points thanks to the proposed decomposition ($(N_x; N_z) = (64; 17)$) while results of original model are still not converged with ($(N_x; N_z) = (256; 65)$). At the same time, a great reduction of dispersion and amplitude errors is noticed. Consequently, in this particular test case, efficiency improves by a reduction between 51 – 79% for a given level of accuracy according to Table 4.1.1, which has to be associated to lower iteration counts with proposed enhancement.

Long-time propagation is a very complex test case for any numerical method as stated in [34] where the original method and a pseudo-spectral High-Order Spectral (HOS) method are compared in terms of accuracy and efficiency. Similar levels of accuracy are achieved with both methods. Particularly, large phase shifts may be noted with low discretization, which is a typical behavior also experienced with other methods, see [35] for instance. Thanks to this test case, accuracy of resolution as well as time integration is assessed. One observes large differences between original method and the enhanced one proposed in this paper. More accurate solution with reduced CPU time is achieved thanks to the explicit description of wavefield with non-linear decomposition as described in Sec. 3.1.

In addition, one may see in the previous table a correlation between phase shift and amplitude error and it appears interesting to study this feature with more details. Indeed, with the enhanced model, if the incident wave field is out of phase with the simulated one, the wave field actually solved will not be in phase with the computed wave field, similarly to what has been observed in Sec. 3.2. Furthermore, it is intended to be associated with larger errors.
In our simulations, in order to keep the phase shift limited, the incident wave phase velocity is chosen to be the numerical one (i.e. obtained from computed solution) instead of the theoretical one from stream function theory. This incident phase velocity is chosen such that incident and computed wave (i.e. sum of incident and scattered components) are in phase. This is ensured by evaluating the phase of the computed solution at each time step and determining the corresponding phase velocity. The tests highlights the importance of selecting an incident wave field which is close to the actual solution. Results are presented in Fig. 12 which presents for the case $N_x = 16$ the evolution during time of phase-shifting $\epsilon_\phi$ and amplitude error $\epsilon_a$. In this figure is compared the original model and the proposed enhanced one with exact and computed phase velocity.

<figure><img src="image_12.png" /><figcaption>Figure 12: Error on the amplitude of the computed wave (right) and phase shift (left) compared to analytical solution. $N_{x\lambda} = 16$, $N_z = 7$, Tolerance is $10^{-6}$ and Courant number $C_r = 0.64$</figcaption></figure>

This figure indicates clearly the enhancement of the non-linear decomposition approach which exhibits lower phase shift errors as well as lower errors on the amplitude of the simulated wave. This enhancement is more
obvious with the phase velocity chosen from the computed solution. Indeed, this enables a limited phase shift and a corresponding low amplitude error on the computed wave field. This figure indicates the strong link between errors on phase shift and amplitude, particularly when non-linear splitting is used. At the same time, it seems interesting to have a look at the algorithmic efficiency with both approaches. This is done in Fig. 13 which presents the evolution during time of the number of iterations needed to solve the linear system.

<figure><img src="image_13.png" /><figcaption>Figure 13: Iteration count during time (GMRES) with non-linear decomposition. $N_{xλ} = 16$, $N_z = 7$, Tolerance is $10^{-6}$ and Courant number $C_r = 0.64$</figcaption></figure>

It is seen that with nonlinear decomposition, at the beginning of the simulation (during ≈ 10 T), the linear system is solved immediately (1 iteration) before an evolution corresponding to the different degrees of complexity of the scattered wave field actually solved by the method. However, the non-linear decomposition appears always more efficient than the original approach on this test case, thanks to information about wave used to setup incident wave
field. It is indeed more accurate (lower phase shift and amplitude errors) and less expensive (in terms of iteration count).

We indicate that in a 3D configuration where wave-structure interactions are studied, the computational domain is usually chosen around the structure of interest (with smallest size possible). Simulations are performed with absorbing conditions of the scattered wave field. Consequently, this wave field and the corresponding errors will not accumulate during the simulation as it happens in this test case (reminding that periodic boundary conditions are used). The induced phase shift is therefore intended to be kept very limited in typical simulations (even over long time propagation), possible errors being propagated and absorbed. Consequently, for typical simulations (equivalently limited duration in periodic condition), high efficiency in the description of incident wave pattern is observed.

### 4.1.2. Convergence rate

We present in this paragraph a study of the convergence rate of the decomposition approach compared to the original OceanWave3D model. The test-case of the reflection of a wave by a vertical wall (see Sec. 3.2) is used here in a linear framework. This allows the definition of an analytical solution to the problem to compare with numerical results. Figure 14 presents the error on the standing wave height at the wall function of the number of points per wavelength used. The error $\epsilon$ is estimated as the relative error between numerical simulation and the analytical exact result: $\epsilon = \frac{|H_{\text{num.}} - H_{\text{exact}}|}{H_{\text{exact}}}$. The expected convergence rate (6th order finite difference scheme is used here) is also shown. The vertical discretization is $N_z = 25$ and numerical parameters are chosen for convergence to be achieved in every simulations (proper
Courant number and absolute tolerance fixed to $10^{-15}$).

<figure><img src="image_14.png" /><figcaption>Figure 14: Relative error $\epsilon$ on standing wave height at wall. Linear simulations compared to analytic result. $N_z = 25$, Courant number $C_r = 0.1$, Tol. = $10^{-13}$</figcaption></figure>

We observe in this figure that both models converge following the expected $6^{th}$ order rate. Convergence is achieved with $N_x = 36$ or $N_x = 24$ points per wavelength on this test-case for decomposition and original approach respectively: accuracy is limited to $10^{-6}$ to $10^{-7}$ because of truncation errors associated with the choice of spatial resolution. Our choice reflects what is often considered a well resolved wave mode (notice that the model is not limited to this accuracy level). For small number of points, the decomposition is more accurate than the original model thanks to the explicit description of the incident wave field. For larger number of points, behavior of both models is similar, including at converged state. The main point shown here is that interesting convergence properties of the original model
are conserved with the new approach.

### 4.1.3. Algorithmic efficiency and computational effort

We detail here an analysis of both algorithmic efficiency (which is independent on implementation) and computational effort (which depends on implementation and hardware) for the non-linear decomposition approach compared to the original OceanWave3D model. Firstly, we present the free propagation case and afterwards the case of the reflection of waves on a vertical wall is studied in details.

Regarding the free propagation case presented in §3.1, the number of iterations needed to solve the linear system is greatly reduced with the decomposition approach. Indeed, as seen in Sec.3.1, for the original model a mean number of iterations of 3.7 was needed while the new approach requires only 1 iteration at each stage (iterative solver is used with a given relative tolerance $10^{-13}$ and an absolute one fixed to $10^{-11}$) during the whole simulation. This comes from the incident wave field which satisfies exactly the system solved. Thus, the residual evaluated in the iterative solver is very low: no scattered waves created during solution process. We remind here that this residual is only due to scattered field in this enhanced model: the system solved (Eqn. 6) is essentially $A\Phi^s = b$ with $A$ a transformed Laplacian and $b$ containing inhomogeneous contributions due to incident field and free surface boundary conditions. This low iteration count is associated with high accuracy resulting in significant reduction of computational effort. These features are conserved with long-time propagation (Sec. 4.1.1) and periodic boundary conditions.
For the case of the reflection of waves on a vertical wall studied in Sec. 3.2 (nonlinear), it appears that the number of iterations needed for inversion is slightly decreased compared to the original method. Comparison of iteration counts and total CPU time for 1000 time steps is presented in Tab. 2 for different number of points, a Courant number $C_r = 0.64$ and a tolerance fixed to $10^{-10}$.

<table><caption>Table 2: Comparison of iteration count and CPU time for the nonlinear reflection on a wall</caption><tbody><tr><th></th><th></th><th>av. iter.</th><th>total CPU time (s)</th></tr><tr><td rowspan="3">N<sub>x</sub> = 96, N<sub>z</sub> = 10</td><td>Original</td><td>12.0</td><td>9.1</td></tr><tr><td>Enhanced</td><td>11.5</td><td>11.6</td></tr><tr><td>Diff.</td><td>-3.7 %</td><td>+27 %</td></tr><tr><th></th><th></th><th>av. iter.</th><th>total CPU time (s)</th></tr><tr><td rowspan="3">N<sub>x</sub> = 192, N<sub>z</sub> = 10</td><td>Original</td><td>9.6</td><td>17.0</td></tr><tr><td>Enhanced</td><td>9.0</td><td>21.6</td></tr><tr><td>Diff.</td><td>-5.6 %</td><td>+27 %</td></tr><tr><th></th><th></th><th>av. iter.</th><th>total CPU time (s)</th></tr><tr><td rowspan="3">N<sub>x</sub> = 96, N<sub>z</sub> = 40</td><td>Original</td><td>12.2</td><td>37.5</td></tr><tr><td>Enhanced</td><td>12.1</td><td>40.0</td></tr><tr><td>Diff.</td><td>-0.4 %</td><td>+6.7 %</td></tr></tbody></table>

In this table, it is seen that the computing time is greater for this test case in OceanWave3D with decomposition. This is due to the computational effort necessary to evaluate incident wave field. However, the most interesting
point to focus on is the average number of iterations needed. A lower iteration
count indicates that when one will consider very large number of points in
2D (or in 3D), the incident terms computations will become negligible and
corresponding CPU time will therefore be reduced with the enhanced model
proposed. This trend is clearly seen in the case $N_x = 96, N_z = 40$ which
exhibits almost same CPU time between both approaches.

Once this is noted, it is important to remind that the considered test case
is probably the most complicated one for the approach of interest. Indeed,
as explained previously, it solves a very complex scattered field due to the
phase shifting between incident wave field and solution of this wave-structure
interaction problem. Consequently, if in this configuration the method ex-
hibits improved numerical features, we are confident that this behaviour can
be extrapolated to any configuration. Furthermore, with typical 3D wave-
structure interactions, it is assumed that this enhancement will be better and
associated with a reduced size of computational domain with this decompo-
sition approach. This will lead to large reduction of computational effort,
reminding that convergence properties of the original model are conserved.

Finally, we present the results of the numerical and algorithmic efficiency
required to solve this problem of waves reflecting off the wall in a specific
configuration. One tries to combine two previous features to define an appro-
priate numerical configuration. These are: i) accurate description of incident
wave field with few points per wavelength with decomposition and ii) com-
parable accuracy between both approaches for the treatment of interaction
between waves and a vertical wall. Therefore, we choose to perform nu-
merical simulations with a refined mesh near the left wall for the enhanced proposed model and a fine regular grid for the original model. The stretching in horizontal direction is achieved thanks to the following definition of grid points : $x_i = L_x \cdot (1 - \frac{\tanh(C_s s_i)}{\tanh(C_s)})$, with $s_i = (N_x - i)ds$, $ds = 1/(N_x - 1)$ and $C_s$ a constant that controls the stretching (chosen as $C_s = 1$ in the following). The comparison in terms of accuracy between both approaches is done in Fig. 15. The problem is solved in a linear framework to have an analytical reference solution.

<figure><img src="image_15.png" /><figcaption>Figure 15: Relative error $\epsilon$ on standing wave height at wall. Linear simulations compared to analytic result. $N_z = 25$, Courant number $C_r = 0.1$, Tol. $= 10^{-13}$</figcaption></figure>

It appears clearly in this figure that the refinement of grid allows a great enhancement of accuracy for a given number of point. Identically, for a given level of accuracy, the number of points needed with the enhanced model is less than for the original model. For this configuration, around half the
number of points is necessary to ensure the similar accuracy between the two approaches. Figure 16 presents an example of meshes with similar accuracy: a stretched grid with 8 points per wavelength with enhanced model and an even grid with 16 points per wavelength for the original model.

<figure><img src="image_16.png" /><figcaption>Figure 16: Example of meshes $N_z = 10$. Top: stretched grid $N_x \frac{\lambda}{L_x} = 8$; bottom: even mesh $N_{x_\lambda} = 16$.</figcaption></figure>

The grid size is similar close to the wall between the two approaches (which ensures the similar accuracy). However, far from the wall, only a few points are needed with the non-linear decomposition (incident wave field explicitly known and scattered waves not of interest in this zone). The chosen set-up for this 2D case is intended to reproduce more precisely the behavior of the non-linear decomposition approach during 3D wave-structure interaction computations. The improvement in terms of computing time is given in Table 3. For this test case rel. and abs. tolerance are fixed to $10^{-6}$ and $10^{-4}$ respectively.
<table><caption>Table 3: Improvement of CPU time with enhanced model (reflection on a wall)</caption><thead><tr><th></th><th></th><th>total CPU time (s)</th><th>av. iter.</th></tr></thead><tbody><tr><td rowspan="3">Case 1, N<sub>z</sub> = 10</td><td>Original N<sub>xλ</sub> = 16</td><td>10.0</td><td>1.16</td></tr><tr><td>Enhanced N<sub>x</sub> λ / L<sub>x</sub> = 8</td><td>4.56</td><td>1.16</td></tr><tr><td>Diff.</td><td>-54 %</td><td>=</td></tr><tr><td></td><td></td><td>total CPU time (s)</td><td>av. iter.</td></tr><tr><td rowspan="3">Case 2, N<sub>z</sub> = 10</td><td>Original N<sub>xλ</sub> = 32</td><td>38.4</td><td>1.15</td></tr><tr><td>Enhanced N<sub>x</sub> λ / L<sub>x</sub> = 16</td><td>17.3</td><td>1.15</td></tr><tr><td>Diff.</td><td>-55 %</td><td>=</td></tr><tr><td></td><td></td><td>total CPU time (s)</td><td>av. iter.</td></tr><tr><td rowspan="3">Case 3, N<sub>z</sub> = 40</td><td>Original N<sub>xλ</sub> = 16</td><td>33.7</td><td>1.16</td></tr><tr><td>Enhanced N<sub>x</sub> λ / L<sub>x</sub> = 8</td><td>15.8</td><td>1.20</td></tr><tr><td>Diff.</td><td>-53 %</td><td>+3 %</td></tr></tbody></table>
This table indicates that the average iteration count is similar between both approaches. Furthermore, thanks to the choice of the numerical set-up which takes advantage of the specific features of the proposed approach, the corresponding computational effort appears greatly improved. The three different set-up indicate that the feature observed is conserved with higher number of points in both directions. In specific 2D configurations (like the one chosen here) or 3D ones in which we study wave-structure interactions, the computational effort will be greatly improved with the proposed enhanced model. This enables to consider the solution of rather complicated problems with this method which appears to be efficient, accurate and robust.

## 4.2. Variable bathymetry

A possible limitation of the non-linear decomposition method is the solution of problems involving variable bathymetry with constant depth incident wave models. In the following, the test case of a wave shoaling on a beach is studied. This is assumed to be the most complicated configuration with this concern. This problem will be solved in a linear framework and compared to the results presented in Bingham & Zhang [2].

### 4.2.1. Set-up of the shoaling wave problem

Following Bingham & Zhang [2], we use a bathymetry defined by

$$h(x) = h_0 - \frac{h_0 - h_1}{2} \left[ 1 + \tanh \left( \frac{\sin \left( \frac{\pi x}{l} \right)}{1 - \left( \frac{2x}{l} \right)^2} \right) \right], \quad -\frac{l}{2} \le x \le \frac{l}{2} \quad (11)$$

with $h_0 = 0.5$, $h_1 = 0.025$ and $l = 12$. This gives a mean slope of $h_x \simeq 0.05$. The corresponding bathymetry profile is given in Fig. 17. Absorbing zones are defined at both sides of the domain ($0 \le x/\lambda \le 2$) and ($14 \le x/\lambda \le 16$).
The waves are propagating from left (i.e. deep water part with a wavelength $\lambda = 1m$ and a corresponding period $T = 0.8018s$) to the right.

<figure><img src="image_17.png" /><figcaption>Figure 17: Bottom profile</figcaption></figure>

### 4.2.2. Results

One expects that this is a difficult case for the new approach, since we define the incident wave field to be a deep water wave with a wavelength $\lambda_0 = 1m$ everywhere in the domain. The scattered wave field will then take into account the presence of the beach and the resulting change in wavelength (at $h_1 = 0.025$, $\lambda_1 \simeq 0.39m$). Figure 18 plots the incident, scattered and total wave field elevation at a given time step once a steady state has been reached (note the different y-axis range for the scattered wave field).

The figure confirms the complexity of the proposed test-case for the decomposition approach. The scattered wave field $\zeta^S$ appears complicated and consequently the solution process is not as efficient as the original one in this case looking at the iteration count (see Tab. 4).
<figure><img src="image_18.png" /><figcaption>Figure 18: Decomposition method: Wave profile of linear wave shoaling up a beach. Top left, incident waves $ζ^I$; Top right, scattered waves $ζ^S$ and Bottom, total wave field $ζ^I + ζ^S$</figcaption></figure>
To ensure the correct solution, the calculations are now compared to both the original approach and the exact solution from energy conservation:

$$ \frac{a_{ref}(x)}{a_0} = \left( \frac{k(x) \left[ 1 + \frac{2k_0h_0}{\sinh(2k_0h_0)} \right]}{k_0 \left[ 1 + \frac{2k(x)h(x)}{\sinh(2k(x)h(x))} \right]} \right)^{\frac{1}{2}} \quad (12) $$

Figure 19 plots the amplitude envelope $a(x)$ of the converged steady-state time signals for the case resolved with 32 points per wavelength (deep water). The analytic solution and the numerical one are in excellent agreement, which ensures the correct solution of the problem. Note that the result obtained with the original model is not given in this figure for clarity since it cannot be distinguished visually from the numerical solution with decomposition.

<figure><img src="image_19.png" /><figcaption>Figure 19: Harmonic analysis of a linear wave shoaling up a beach</figcaption></figure>

We now consider the computational effort required to solve the problem by comparing the number of iterations per Runge-Kutta stage. Results are
reported in Tab. 4. In this table is also reported the relative error $\epsilon$ between numerical solutions $a_{num}$ and analytical results $a_{ref}$ defined as:

$$ \epsilon = \frac{\int_{2}^{14} |a_{num}(x) - a_{ref}(x)| dx}{\int_{2}^{14} |a_{ref}(x)| dx} $$

The numerical simulations are performed with the exact same numerical parameters: absorption zone, Courant number $C_r = 1$, Tolerance (rel. $10^{-9}$ and abs. $10^{-7}$).

<table><caption>Table 4: Evolution of iteration count with number of points per wavelength $N_{xλ}$. Rel. tolerance $10^{-9}$ and abs. tolerance $10^{-7}$</caption><thead><tr><th rowspan="2"></th><th colspan="2">N<sub>xλ</sub> = 8</th><th colspan="2">N<sub>xλ</sub> = 16</th><th colspan="2">N<sub>xλ</sub> = 32</th><th colspan="2">N<sub>xλ</sub> = 64</th></tr><tr><th>iter</th><th>ε</th><th>iter</th><th>ε</th><th>iter</th><th>ε</th><th>iter</th><th>ε</th></tr></thead><tbody><tr><th>Original</th><td>6.4</td><td>1.1 10<sup>-1</sup></td><td>4.3</td><td>3.6 10<sup>-3</sup></td><td>3.3</td><td>3.6 10<sup>-4</sup></td><td>2.7</td><td>7.8 10<sup>-4</sup></td></tr><tr><th>Enhanced</th><td>5.8</td><td>3.6 10<sup>-2</sup></td><td>4.7</td><td>1.8 10<sup>-3</sup></td><td>3.8</td><td>1.7 10<sup>-4</sup></td><td>2.8</td><td>1.8 10<sup>-4</sup></td></tr><tr><th>Diff.</th><td>-9.4%</td><td>-67%</td><td>+9.3%</td><td>-50%</td><td>+15%</td><td>-52%</td><td>+3.7%</td><td>-77%</td></tr></tbody></table>

We observe that the required number of iterations is slightly higher with the new approach, reminding that the proposed test-case is severe for the method. At the same time, the accuracy appears to be greatly enhanced with the decomposition approach (notice that the error is larger for $N_{xλ} = 64$ with both approaches due to some reflections in the absorbing zone). Thus, for a given accuracy, computational effort is (at least) equivalent between both approaches. Consequently, even in this severe test-case, the overall computational effort is enhanced with the nonlinear decomposition. It has to be noted that this specific configuration may be solved more efficiently if one uses a dedicated incident wave model able to simulate the influence
of bathymetry. This way, similarly to Sec. 3.1, the incident wave model satisfies all boundary conditions. Then, the previous issues are now obsolete, the enhanced solver just check that the incident wave field satisfies the non-linear system Eqn. (6) to the level of truncation errors. This solution will be preferred when treating wave-structure interactions in domains with variable depth. However, if such model is not available, this test-case indicates that the non-linear decomposition approach is actually able to take into account the influence of the bathymetry accurately and efficiently.

## 5. Conclusion

In this paper the development of an efficient numerical model for solving wave-structure interaction problems in the coastal and ocean environment is described. The model is based on the high-order finite differences solver OceanWave3D (Bingham & Zhang [2] and Engsig Karup et al. [3, 4]). This model solves the Laplace problem over the whole fluid domain and has been shown to be efficient and highly accurate. To enhance the solution process, we adopt a non-linear decomposition of the wave-structure interaction solution into an incident and a scattered wave field, following [17]. The incident part represents the incoming wave which is known explicitly while the scattered part takes into account the presence of the structure inside the given incident wave field. The definition and/or generation of waves (including complex ones) will be easier (especially in 3D configurations) and strong gains in computing time are expected. Indeed, the proposed method allows the use of refined grid around structures for scattering problems. This possibly reduces problem size and consequently enhances efficiency. Different
validations in 2D have been performed, which assessed the robustness and accuracy of the proposed decomposition approach. Particularly, the extreme case of non-linear waves reflecting on a vertical wall has been studied in details with favorable comparison to the original OceanWave3D model. The steep standing waves appearing close to the wall are analysed and especially the steepest stable one which is successfully compared to experiments and to a dedicated numerical model. Finally, the efficiency of the non-linear decomposition is studied. The accuracy of the model to propagate wave fields during long time has been assessed. Then, a study of convergence rate indicates that the decomposition approach conserves the interesting convergence properties of the original model. Those features are used to elaborate a test-case which confirms the effective performance benefits one may obtain with the proposed method in 2D and 3D configurations. The last section indicates that for the depth-varying computational domains, the use of the proposed model leads also to enhanced computational effort even if a specific incident wave model may be used instead of the classical constant-depth one for optimal efficiency. The resulting method proposed in this paper is consequently robust, accurate and efficient.

As a conclusion, the different studied 2D test cases presented the main advantages of the proposed approach: i) improved computational efficiency in solve step, ii) reduction of problem size for localized scattering, iii) conservation of convergence properties (w.r.t original OceanWave3D model) and iv) possibility to treat depth varying test cases. The chosen test cases are supposed to be the most severe for the decomposition approach and their detailed analysis proposed in this paper exhibit enhanced computational effort.
Consequently, significant gains can be achieved on, a priori, any configurations. This method will be particularly efficient in 3D wave-structure interactions where scattered field is localized around the structure. The extension to 3D bottom-mounted geometries is straightforward with promising results obtained in Ducrozet et al. [23]. A follow up paper (part B) is intended to be published which highlights the efficiency of the methodology. The next step will be to take into account more complex geometries (in order to treat floating structures for instance). The use of multi-block, boundary-fitted solutions seems adequate with encouraging preliminary results. On the other hand, other incident wave models can be used to take into account complex wave fields [24, 36] (e.g. irregular 3D sea states and extreme waves).

## Acknowledgments

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

[1] A. Di Mascio, M. Landrini, F. Lalli, U. Bulgarelli, Three dimensional non linear diffraction around fixed structures, in: Proceedings of the 20th Symposium on Naval Hydrodynamics.

[2] H. B. Bingham, H. Zhang, On the accuracy of finite-difference solutions for nonlinear water waves, J. Eng. Math. 58 (2007) 211–228.
[3] A. P. Engsig-Karup, H. B. Bingham, O. Lindberg, An efficient flexible-order model for 3d nonlinear water waves, J. Comp. Phys. 228 (2009) 2100–2118.

[4] A. P. Engsig-Karup, M. G. Madsen, S. L. Glimberg, A massively parallel gpu-accelerated model for analysis of fully nonlinear free surface waves, International Journal for Numerical Methods in Fluids 70 (2012) 20–36.

[5] Š. Malenica, B. Molin, Third harmonic wave diffraction by a vertical cylinder, J. Fluid Mech. 302 (1995) 203–229.

[6] M. Celebi, M. Kim, R. Beck, Fully nonlinear three-dimensional numerical wave tank simulation, Journal of Ship Research 42 (1998) 33–45.

[7] G. Weymouth, R. Wilson, F. Stern, RANS computational fluid dynamics predictions of pitch and heave ship motions in head seas, J. Ship Research 49 (2005) 80–97.

[8] M. S. Longuet-Higgins, E. D. Cokelet, The deformation of steep surface waves on water. i. a numerical method of computation, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 350 (1976) 1–26.

[9] W. Bai, R. Eatock Taylor, Fully nonlinear simulation of wave interaction with fixed and floating flared structures, Ocean Engineering 36 (2009) 223 – 236.

[10] M. Christou, C. H. Hague, C. Swan, The reflection of nonlinear irregular surface water waves, Engineering Analysis with Boundary Elements 33 (2009) 644 – 653.
[11] C. Fochesato, S. Grilli, F. Dias, Numerical modeling of extreme rogue waves generated by directional energy focusing, Wave Mation 44 (2007) 395–416.

[12] G. Wu, R. Eatock Taylor, Finite element analysis of two dimensional non-linear transient water waves, Appl. Ocean Res. 16 (1994) 363–372.

[13] Q. Ma, G. Wu, R. Eatock Taylor, Finite element simulation of fully non-linear interaction between vertical cylinders and steep waves - part 1: methodology and numerical procedure, Int. J. Num. Methods Fluids 36 (2001) 265–285.

[14] B.Li, C. Fleming, A three dimensional multigrid model for fully nonlinear water waves, Coastal Engng. 30 (1997) 235–258.

[15] S. L. Glimberg, A. P. Engsig-Karup, M. G. Madsen, A fast gpu-accelerated mixed-precision strategy for fully nonlinear water wave computations, in: Proceedings of European Numerical Mathematics and Advanced Applications (ENUMATH).

[16] P. Ferrant, Simulation of strongly-nonlinear wave generation and wave-body interaction using a 3D MEL model, in: Proc. 21st ONR Symp on Naval Hydrodynamics, Trondheim, Norvège, pp. 93–109.

[17] P. Ferrant, Fully non-linear interactions of long-crested wave packets with a three-dimensional body, in: Proc. 22nd ONR Symposium on Naval Hydrodynamics, Washington, USA, pp. 403–415.

[18] P. Ferrant, L. Gentaz, B. Alessandrini, D. Le Touzé, A potential /
RANSE approach for regular water-wave diffraction about 2d structures,
Ship Tech. Res. 50 (2003).

[19] R. Luquet, B. Alessandrini, P. Ferrant, L. Gentaz, Simulation of the viscous flow past a ship in waves using the swense approach, in: Proc. 24th ONR Symposium on Naval Hydrodynamics, St John's, Newfoundland.

[20] R. Luquet, G. Ducrozet, L. Gentaz, P. Ferrant, B. Alessandrini, Applications of the swense method to seakeeping simulations in irregular waves, in: Proc of the 9th Int. Conf. on Num. Ship Hydro., Ann Arbor, Michigan.

[21] V. Zakharov, Stability of periodic waves of finite amplitude on the surface of a deep fluid, J. Appl. Mech. Tech. Phys. (1968) 190-194.

[22] J. D. Fenton, The numerical solution of the steady water wave problem, Comp. & Geosc. 14 (1988) 357-368.

[23] G. Ducrozet, H. B. Bingham, A. P. Engsig-Karup, P. Ferrant, High-order finite difference solution for 3d nonlinear wave-structure interaction, Journal of Hydrodynamics, Ser. B 22 (2010) 225-230.

[24] G. Ducrozet, F. Bonnefoy, D. Le Touzé, P. Ferrant, A modified high-order spectral method for wavemaker modeling in a numerical wave tank, European Journal of Mechanics - B/Fluids 34 (2012) 19-34.

[25] A. Savitsky, M. Golay, Smoothing and differentiation of data by simplified least squares procedures, Anal. Chem. 36 (1964) 1627-1639.
[26] J. E. Kuo, H. Wang, S. Pickup, Multidimensional least-squares smooth-
ing using orthogonal polynomials, Anal. Chem. 63 (1991) 630–635.

[27] J. Berland, C. Bogey, O. Marsden, C. Bailly, High-order, low dispersive and low dissipative explicit schemes for multiple-scale and boundary problems, J. Comp. Phys. 224 (2007) 637–662.

[28] C. H. Su, R. M. Mirie, On head-on collisions between two solitary waves,
J. Fluid Mech. 98 (1980) 509–525.

[29] J. D. Fenton, M. M. Rienecker, A fourier method for solving nonlinear water wave problems: application to solitary-wave interactions, J. Fluid Mech. 118 (1982) 411–443.

[30] G. Mercer, A. Roberts, Standing waves in deep water: Their stability and extreme form, Phys. Fluids A 4 (1992) 259–269.

[31] C. Tsai, D. Jeng, Numerical fourier solutions of standing waves in finite water depth, Appl. Ocean Res. 16 (1994) 185–193.

[32] M. Okamura, Standing gravity waves of large amplitude in deep water,
Wave Motion 37 (2003) 173–182.

[33] G. I. Taylor, An experimental study of standing waves, Proc. Roy. Soc.
A 218 (1953) 44–59.

[34] G. Ducrozet, H. B. Bingham, A. P. Engsig-Karup, F. Bonnefoy, P. Fer-
rant, A comparative study of two fast nonlinear free-surface water wave
models, International Journal for Numerical Methods in Fluids 69 (2012)
1818–1834. Doi: 10.1002/fld.2672.
[35] D. Fructus, D. Clamond, J. Grue, O. Kristiansen, An efficient model for three-dimensional surface wave simulations. Part I: Free space problems, J. Comp. Phys. 205 (2005) 665–685.

[36] G. Ducrozet, F. Bonnefoy, D. Le Touzé, P. Ferrant, 3-D HOS simulations of extreme waves in open seas, Nat. Hazards Earth Syst. Sci. 7 (2007) 109–122.