M Gouin

G Ducrozet

P Ferrant

Development and validation of a non-linear spectral model for water waves over variable depth

Keywords: Hydrodynamics, Non-linear water waves, High-Order Spectral method, Variable bathymetry

In the present paper two numerical schemes for propagating waves over a variable bathymetry in an existing High-Order Spectral (HOS) model are introduced. The first scheme was first developed by Liu and Yue (1998), and the second one is an improved scheme which consider two independent orders of non-linearity: one for the bottom and one for the free-surface elevation. We investigate the numerical properties (accuracy, convergence rate, efficiency) of both schemes with respect to the numerical parameters on a simple configuration. To validate the proposed schemes, we first consider Bragg reflection from a sinusoidal bottom patch -as an example of a small bottom variation around a mean water depth. The second validation case focuses on a larger bottom variation with the study of the shoaling of linear waves. Finally, an application is performed to demonstrate the applicability of the proposed model to non-linear cases where the bottom variation is important. In this concern, the very well documented test case of a 2D underwater bar is studied in details. Comparisons are provided with both experimental and numerical results.

Introduction

The accurate modelling of surface gravity waves over non-negligible bottom topography is of major interest in ocean engineering, especially in the field of marine renewable systems. These marine structures are intended to be deployed in limited water depth, where the effect of variable bathymetry on local wave conditions is very significant. Indeed, when entering shallow water zones, waves are strongly affected by the bottom through shoaling, refraction, diffraction, reflection and the resulting variations in local wave speed. Thus, the accurate description of the wave field over variable depth is a prerequisite for the accurate prediction of wave loads acting on structures in coastal zones.

For this purpose, a wide variety of non-linear flow models have been developed during the last decades. Some of them are based on the solution of the Reynolds Averaged Navier-Stokes equations, such as the CFD models presented in [START_REF] Westphalen | Investigation of wave-structure interaction using state of the art CFD techniques[END_REF], but the computational effort with these models remains prohibitive. Thus, most of the non-linear flow models for wave propagation were developed in the framework of the potential flow theory, considering that the propagation in the ocean is mostly irrotational and inviscid (neglecting wave breaking at the sea surface and dissipation due to bottom friction).

The Boundary Element Model (BEM) (see for instance [START_REF] Grilli | A fully non-linear model for threedimensional overturning waves over an arbitrary bottom[END_REF] is one of the methods used to represent wave propagation over non-uniform depth in wide domains. The problem is solved on the boundaries, allowing the reduction of the problem size. Moreover, it is a very flexible method as it can account for a variable bottom and very complex geometries, including structures. Nevertheless, it requires the inversion of full matrices, reducing the efficiency of the method. Some recent developments using the Fast Multiple Acceleration technique (see [START_REF] Korsmeyer | Multipoleaccelerated preconditioned iterative methods for three-dimensional potential problems[END_REF]; [START_REF] Fochesato | A fast method for nonlinear three-dimensional free-surface waves[END_REF] intend to improve the efficiency of the BEM model. Finite difference methods (see [START_REF] Engsig-Karup | An efficient flexible-order model for 3D nonlinear water waves[END_REF]) are also flexible in terms of geometries, but as any volume-type method it requires a high number of nodes to represent the whole computational domain. However, this formulation leads to sparse matrices allowing the use of advanced numerical procedures and resulting in a good efficiency. Other methods for modelling wave propagation over variable depth can also be found in the literature such as Finite Element Methods [START_REF] Wu | Finite element analysis of two-dimensional nonlinear transient water waves[END_REF]; [START_REF] Ma | Quasi ALE finite element method for nonlinear water waves[END_REF]) or a recent fully dispersive coupled-mode model described by [START_REF] Belibassakis | Directional wave spectrum transformation in the presence of strong depth and current inhomogeneities by means of coupled-mode model[END_REF]. Boussinesq methods were initially developed for small relative water depths in the framework of weakly non-linear waves, but the last developments of high order versions of Boussinesq approximations developed by [START_REF] Bingham | Nonlinear wave interaction with bottom-mounted structures by a high-order Boussinesq method[END_REF] and [START_REF] Madsen | A Boussinesq-type method for fully nonlinear waves interacting with a rapidly varying bathymetry[END_REF] allow to account for larger water depths at little extra computational cost which remains quite attractive.

Interesting properties of spectral methods in terms of convergence have led to the development of numerous models. We can first cite the Direct Method, introduced by [START_REF] Fenton | A Fourier method for solving nonlinear waterwave problems: application to solitary-wave interactions[END_REF], which solves the problem on the free-surface (free-surface method) allowing a reduction of one dimension. Nevertheless, the required inversion of a fully populated matrix makes it not very efficient. The pseudo-spectral σ-transform model introduced by [START_REF] Chern | A pseudospectral σtransformation model of 2-D nonlinear waves[END_REF] allows the modelling of more complex geometries but still with a high computational time due to the necessity to discretize the whole fluid domain.

The pseudo-integral/spectral method improved by [START_REF] Fructus | An efficient model for three-dimensional surface wave simulations[END_REF] uses a pseudo-spectral solution added to an integral solution to cope with steeper cases and variable bathymetries. The DNO (Dirichlet to Neumann Operator) method was initiated by [START_REF] Craig | Numerical simulation of gravity waves[END_REF] for a flat bottom. This method was extended to a variable bottom by [START_REF] Guyenne | Numerical simulation of solitary waves on plane slopes[END_REF] and [START_REF] Craig | Hamiltonian long-wave expansions for water waves over a rough bottom[END_REF] by introducing another operator depending only on the variation of the bottom. It was then improved numerically by [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF].

In the present paper, we use the High-Order Spectral (HOS) method. This non-linear potential method has been initially developed by [START_REF] West | A new numerical method for the surface hydrodynamics[END_REF] and [START_REF] Dommermuth | A high-order spectral method for the study of non-linear gravity waves[END_REF] for a flat bottom, and extensively validated for different configurations in the LHEEA Lab., from regular waves up to irregular multidirectional wavefields (see [START_REF] Ducrozet | 3-D HOS simulations of extreme waves in open seas[END_REF]Ducrozet et al. ( , 2012b))). This model, named HOS-ocean is available as an open-source version1 . The HOS formalism presents expansions identical to the DNO method, as demonstrated in [START_REF] Schäffer | Comparison of Dirichlet-Neumann operator expansions for nonlinear surface gravity waves[END_REF], and several advantages. First of all, the problem is formulated on the free surface, allowing a reduction of one dimension when solving the problem.

Moreover, it allows the solution of the problem with the fully non-linear free surface boundary conditions, it shows excellent convergence properties and it has a low computational cost. Therefore this method is very efficient and accurate, but is initially limited to simple geometries in both horizontal and vertical directions. Non-linear potential flow models cited above consider a varying bottom but few of them were based on the HOS scheme. As demonstrated in Ducrozet et al. (2012a), the HOS model appears more efficient than the most advanced potential flow solvers, when compared for wave propagation on uniform depth.

Thus it appears interesting to extend such a model to a variable bathymetry in order to broaden its application range, while keeping the numerical efficiency. [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF] provided one study with the HOS method which takes into account a variable bathymetry, and presented one simulation case but limited to a small variation of the bottom. The first scheme used in the present paper is based on their work. The second scheme is an improvement of the first one by considering two different orders of non-linearity: one for the bottom and one for the free-surface. This dissociation of the orders of non-linearity was presented with the DNO method in [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF], and will be adapted here to the HOS formalism.

The two HOS methods presented (called original method and improved method) are explained in detail in the present paper. Some of the work with the original method presented hereafter has already been partially introduced in [START_REF] Gouin | Development and validation of a highly nonlinear model for wave propagation over a variable bathymetry[END_REF]. For documenting the accuracy of both schemes, a test case with a constant variation of the bottom is computed to check the convergence on a simple geometry with a non-linear reference solution. It is a highly demanding test case because a constant but wrong depth is imposed in the whole domain. Then both methods are applied to 2D monochromatic cases.

Bragg reflection from a sinusoidal bottom patch will first be considered, as an example of a small bottom variation around a mean water depth. To validate our models with larger bottom variation, the shoaling of linear waves is studied in the second case. Finally, two well-documented application cases (both numerically and experimentally) are provided. We validate and compare the two methods for realistic and large bottom variations. The first case considers the transformation of a non-linear, monochromatic wave as it travels up and over a submerged bar with a mild slope. This validation case has often been used as a discriminating test case for non-linear models of surface waves propagation over a variable bottom because it requires the accurate propagation of waves in both deep and shallow water. A comparison to the experimental data and to other numerical results is provided, as well as a comparison of both methods.

The second application case considers the same experimental set-up but with a steeper slope, and proves the ability of our models to represent cases with strong variations of the bathymetry and large bottom gradients.

Methods and algorithms

Hypothesis and formulation of the problem

A 2D rectangular fluid domain and a Cartesian coordinate system with the origin O located at one corner of the domain are considered. The z axis is vertical and oriented upwards, with the level z = 0 corresponding to the mean water level. z = η (x, t) represents the free surface elevation, h the total water depth, h 0 the mean depth and β (x) the bottom variation, such as -h (x) = -h 0 +β (x).

Thus, the domain considered is: -h 0 + β (x) ≤ z ≤ η (x) (see Fig. 1). Following [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF], both kinematic and dynamic non-linear freesurface boundary conditions (FSBC) are written in terms of surface quantities η and φ ( φ (x, t) = φ (x, z = η, t)), and expressed at the exact free-surface position

z = η: ∂η ∂t = 1 + ∂η ∂x 2 ∂φ ∂z - ∂ φ ∂x . ∂η ∂x (2) ∂ φ ∂t = -gη - 1 2 ∂ φ ∂x 2 + 1 2 1 + ∂η ∂x 2 ∂φ ∂z 2 (3)
To account for the time evolution of the quantities of interest η and φ one only needs to evaluate the vertical velocity at the free surface W (x, t) = ∂φ ∂z (x, z = η (x, t) , t). This is the purpose of the HOS method described in 2.2.

125

The HOS method was initially developed for a flat bottom but here for a non-constant bathymetry the bottom boundary condition reads :

∂φ ∂x ∂β ∂x - ∂φ ∂z = 0 on z = -h 0 + β (x) (4)
To account for a varying bathymetry, [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF] introduces an additional potential. Thus, the total potential φ tot solution of the problem (which was named φ before) is expressed as:

φ tot = φ h0 + φ β (5)
φ h0 satisfies a Neumann condition on z = -h 0 , therefore φ h0 is solution of the problem at constant depth h 0 :

∂φ h0 ∂z (x, z = -h 0 , t) = 0 on z = -h 0 . (6)
φ β allows the definition of the correct bottom boundary condition (Eq.4)

and satisfies a Dirichlet condition on z = 0:

φ β (x, z = 0, t) = 0 on z = 0. ( 7
)
In 2D, the potentials are expanded on basis functions taking into account previous boundary conditions and the periodicity of the domain:

φ h0 (x, z, t) = j A j (t) cosh (k j (z + h 0 )) cosh (k j h 0 ) e ikj x (8) φ β (x, z, t) = j B j (t) sinh (k j z) cosh (k j h 0 ) e ikj x (9)
with k j = j 2π Lx and A j and B j the modal amplitudes of φ h0 and φ β respectively.

High-Order-Spectral Method

The HOS model is a pseudo-spectral model initially developed in [START_REF] West | A new numerical method for the surface hydrodynamics[END_REF]; [START_REF] Dommermuth | A high-order spectral method for the study of non-linear gravity waves[END_REF]. The potential is expressed as a truncated power series of components φ (m)

tot for m = 1 to M (M is the order of nonlinearity of the HOS method). Then, the potential evaluated at the free surface is expanded in a Taylor series with respect to the mean water level z = 0 (Eq.10).

Combining these two expansions gives a triangular set of Dirichlet problems for the components that can be solved by means of a spectral method (allowing the use of Fast Fourier Transforms (FFTs) for efficient computations).

φtot = φ tot (x, z = η, t) = M m=1 M -m n=0 η n n! ∂ n φ (m) tot ∂z n (x, z = 0, t) (10)
∂φ tot ∂z (x, z = η, t) = W (x, t) = M m=1 m-1 k=0 η k k! ∂ k+1 φ (m-k) tot ∂z k+1 (x, 0, t) (11)
Equations 11 can be reduced to Eq.(5.7) cited in [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF] by an appropriate change of variables.

One more equation (the bottom condition) is needed to find the modal amplitudes of order m : A (t) are computed (see Fig. 2), by evaluating the vertical derivative of the potential and using the same kind of double expansion, the vertical velocity W at the free surface can be obtained. Equations (2,3) provide the time derivatives of the unknowns η and φ tot (as presented in Fig. 3), which are used in a time-marching classical fourth-order Runge-Kutta scheme with an adaptative time step (see [START_REF] Fructus | An efficient model for three-dimensional surface wave simulations[END_REF] for the details of the procedure).

Thus, the free surface conditions (Eqs. 2 and 3) are marched in time in a similar manner as for the initial HOS procedure (see [START_REF] Ducrozet | 3-D HOS simulations of extreme waves in open seas[END_REF]).

Note that a careful dealiasing is performed as explained in [START_REF] Bonnefoy | Time-domain simulation of nonlinear water waves using spectral methods[END_REF] to ensure the method's convergence and accuracy, even for high steepnesses.

Indeed, for high non-linearities, a full or total dealiasing allows us to reduce the errors observed by [START_REF] Dommermuth | A high-order spectral method for the study of non-linear gravity waves[END_REF]. The two following sections will present the addition due to the consideration of a non-uniform bottom.

Original method

The way of using the bottom condition to find the modal amplitudes B The bottom condition (Eq.4) truncated at the order M reads:

M m=1 ∂φ (m) tot ∂x ∂β ∂x - M m=1 ∂φ (m) tot ∂z = 0 on z = -h 0 + β (x) (12)
By assuming that β h0 << 1 we can write a Taylor expansion with respect to the mean depth z = -h 0 :

M m=1 ∂φ (m) tot ∂x ∂β ∂x -h0 + β ∂ ∂z ∂φ (m) tot ∂x ∂β ∂x -h0 + ... + β M M ! ∂ M ∂z M ∂φ (m) tot ∂x ∂β ∂x -h0 = M m=1 ∂φ (m) tot ∂z -h0 + β ∂ 2 φ (m) tot ∂z 2 -h0 + ... + β M M ! ∂ M +1 φ (m) tot ∂z M +1 -h0 (13)
Since we follow the potential theory (Eq.1) we have in two dimensions:

∂ 2 φtot ∂z 2 = -∂ 2 φtot
∂x 2 , and we obtain:

M m=1 ∂φ (m) tot ∂x ∂β ∂x -h0 + ... + β M (M )! ∂ M ∂z M ∂φ (m) tot ∂x ∂β ∂x -h0 = M m=1 ∂φ (m) tot ∂z -h0 -β ∂ 2 φ (m) tot ∂x 2 -h0 -... - β M M ! ∂ M -1 ∂z M -1 ∂ 2 φ (m) tot ∂x 2 -h0 (14)
We also assume O (β) ≡ O ∂β ∂x ≡ O (η) and we only keep terms of order η (m) :

m-1 l=1 β l-1 (l -1)! ∂ l-1 ∂z l-1 ∂φ (m-l) tot ∂x ∂β ∂x -h0 = ∂φ (m) tot ∂z -h0 - m-1 l=1 β l l! ∂ l-1 ∂z l-1 ∂ 2 φ (m-l) tot ∂x 2 -h0 (15)
Thus we find the equations presented in [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF]:

m = 1 : ∂φ (1) tot ∂z (x, -h 0 , t) = 0 m = 2, ..M : ∂φ (m) tot ∂z (x, -h 0 , t) = m-1 l=1 ∂ ∂x β l l! ∂ l-1 ∂z l-1 ∂φ tot ∂x (m-l) z=-h0 (16)
Perturbation expansions are truncated at order M . This development (Eq.( 16)

with Eqs.(8, 9)) allows the computation of the modal amplitudes B j1 (t). Modal amplitudes are independent of the position x, thus FFTs can still be used, preserving the numerical efficiency of the original HOS scheme.

Improved method

As described in [START_REF] Schäffer | Comparison of Dirichlet-Neumann operator expansions for nonlinear surface gravity waves[END_REF], the HOS method exhibits a formalism similar to the DNO method. The improved method is based on the extension of the DNO scheme to a variable bottom proposed by [START_REF] Craig | Hamiltonian long-wave expansions for water waves over a rough bottom[END_REF] and [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF]. Starting from Eq.( 12) and assuming that β h0 << 1 we can write the same Taylor expansion of the bottom boundary condition as before (Eq.14) but at a different order of truncation M b :

M m=1   ∂φ (m) tot ∂x ∂β ∂x -h 0 + β ∂ ∂z ∂φ (m) tot ∂x ∂β ∂x -h 0 + ... + β M b -1 (M b -1)! ∂ M b -1 ∂z M b -1 ∂φ (m) tot ∂x ∂β ∂x -h 0   = M m=1   ∂φ (m) tot ∂z -h 0 -β ∂ 2 φ (m) tot ∂x 2 -h 0 -... - β M b (M b )! ∂ M b -1 ∂z M b -1 ∂ 2 φ (m) tot ∂x 2 -h 0   (17)
In the left-hand side of Eq.17, the expansion is made up to order M b -1 155 because we assume that O (β) ≡ O ∂β ∂x . Thus the left-hand side is of order M b , and is consistent to the right-hand side.

Eq.17 can be written as:

M m=1 ∂φ (m) tot ∂z -h0 = M m =1 M b l=1    ∂ ∂x   β l l! ∂ l-1 ∂z l-1   ∂φ (m ) tot ∂x     -h0    (18)
We adopt the hypothesis that each potential φ

φ (m) β = M b l=1 φ (m,l) β (19)
and the corresponding modal amplitudes are given by:

B (m) j = M b l=1 B (m,l) j (20)
Thus, at each order m, introducing φ

(m) tot = φ (m) h0 + φ (m) β
and keeping only the terms of maximum order m + M b (assuming as for the original method that

O (β) ≡ O ∂β ∂x ≡ O (η)) we get: ∂φ (m) β ∂z -h0 = M b l=1 ∂φ (m,l) β ∂z -h0 = M b l=1 ∂ ∂x β l l! ∂ l-1 ∂z l-1 ∂φ (m) h0 ∂x -h0 + M b l=1 l-1 p=1 ∂ ∂x β p p! ∂ p-1 ∂z p-1 ∂φ (m,l-p) β ∂x -h0 (21)
We can then derive an iterative formula to calculate each B (m,l) j at a given 160 order m. This expression can be found in the DNO formalism in [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF]. Here we use the same formalism to compare our formulas with the ones given in their paper.

For l odd:

kB (m,l) = -k β l l! sech (kh 0 ) k l-1 k -k   l-1 p=2,even β p p! k p-2 kkB (m,l-p) - l-2 p=1,odd β p p! tanh (kh 0 ) k p-2 kkF l-p   + ik ∂β ∂x β l-1 (l -1)! sech (kh 0 ) k l-2 k + ik ∂β ∂x   l-1 p=2,even β p-1 (p -1)! k p-3 kkB (m,l-p) - l-2 p=1,odd β p-1 (p -1)! tanh (kh 0 ) k p-3 kkB (m,l-p)   (22)
And for l even:

kB (m,l) = -k   l-2 p=2,even β p p! k p-2 kkB (m,l-p) - l-1 p=1,odd β p p! tanh (kh 0 ) k p-2 kkB (m,l-p)   + ik ∂β ∂x   l-2 p=2,even β p-1 (p -1)! k p-3 kkB (m,l-p) - l-1 p=1,odd β p-1 (p -1)! tanh (kh 0 ) k p-3 kkB (m,l-p)   (23)
If we compare Eq.( 22,23) to the developments obtained by [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF], one can notice that we obtain an additional part due to the gradient of the bottom (this part is missing in [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF]).

Others terms are identical. Previous equations (22 and 23) should then be used in the DNO method with bottom variation.

In comparison to the original method, we can compute the modal amplitudes

B (m) j
(t) due to the bottom even for m = 1 (with the original method we have to begin at m = 2 for the computation of the

B (m) j (t)).
Moreover, with this improved method, the orders of non-linearity on the free-surface and on the bottom can be selected independently. This gives a better flexibility to the model and should be more convenient to compute cases with different requirements in bottom variation and free-surface non-linearities.

If we choose M b = M with the improved method, the results are a bit different to those obtained with the original method because we do not calculate the same terms at each sub-step i = 1, ..., M . Section 3.1.3 shows the differences on the error made on the vertical velocity with the two methods.

Lateral boundary conditions and relaxation zones

With the present scheme periodic boundary conditions are required in the xdirection. However, practical applications often requires non-periodic boundary conditions (see e.g. the test cases presented). A simple way is to increase the size of the computational domain but it is not suitable for long-time simulations.

Here we proposed to use relaxation zones.

Wave generation and absorption is achieved by relaxing the numerical solution towards a specific solution over regions near the boundary domain. In 2D, these boundary regions are located at the right and left of the domain.

Relaxation zones are created by defining a relaxation coefficient 0 ≤ C r (x) ≤

1 and an exact imposed solution (η imp , φimp ) as described in [START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF] and [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF]. In our cases, the imposed solution will be either the stream solution of [START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF] (for non-linear regular wave definition) or the zero solution (for complete absorption).

After each stage of the time integration in the Runge-Kutta scheme, the solution in the relaxation zone is redefined by:

η (x, t) = (1 -C r (x)) .η imp + C r (x).η (x, t) (24) φ (x, t) = (1 -C r (x)) . φimp + C r (x). φ (x, t) (25)
To avoid reflections, C r should be smooth. A polynomial of order 3 has been used to define it on zones of several wavelengths.

Moreover, when calculating the derivatives of the bottom, FFTs are used, so the bottom needs to be periodic. When computing a non-periodic 2D case, the bathymetry is thus artificially modify at the end of the domain (in the direction of propagation) to match the depth at the beginning of the domain (the same needs to be done in the transverse direction for 3D cases). This does not affect the results as the modification of the bathymetry will take place in the relaxation zones.

Regular wave propagation

A non-linear regular wave propagating over a flat bottom is considered.

The initialisation of η and φ is made by a non-linear regular wave solution of [START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF] built at water depth h. The potential φ h0 is defined at the mean depth h 0 (which differs from h) and the potential φ β must modify and correct the water depth β (β is a constant) which is introduced artificially (cf Fig. 4). The goal is to rebuild the potential φ tot at the total water depth -h (x) = -h 0 + β (x). This validation case is highly demanding because an inadequate depth is imposed in the whole domain. Thus if this test case shows accurate results, other cases will appear less demanding and should also work.

Error on the vertical velocity.

In order to advance in time the unknowns on the free surface η and φ (cf Eqs 2 and 3), we first need to evaluate the vertical velocity W (x, t) = ∂φ ∂z (x, z = η (x, t) , t) associated to φ tot . This test case proposes to study the error w = W -W RF W RF made on the evaluation of the vertical velocity once η and φ are known (W RF denotes the reference vertical velocity from [START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF]). The level of error correctly quantifies the accuracy of the model in the propagation of such waves, as demonstrated in [START_REF] Bonnefoy | Time-domain simulation of nonlinear water waves using spectral methods[END_REF].

The regular wave is described by its steepness ka and the relative water depth kh. The influence of the following parameters is investigated : ka (from 0.01 to 0.2), kh (from 0.5 to 10), the bottom variation β h0 , the number of nodes N and the HOS orders M and M b . Results will be compared to those obtained in the case of a finite constant depth ( β h0 = 0%) which was extensively validated in [START_REF] Bonnefoy | Time-domain simulation of nonlinear water waves using spectral methods[END_REF]. With the original HOS model on flat bottom, the error achieved on the vertical velocity presents an exponential convergence in function of N and M .

This convergence test is used to characterize both methods and has already been presented in [START_REF] Gouin | Development of a highly nonlinear model for wave propagation over a variable bathymetry[END_REF] for the original method. More details and a comparison with the improved method are given hereafter. It allows to assess the convergence of both schemes on the reconstruction of the vertical velocity with a wide variety of wave conditions and non-negligible bottom variations.

Original method

We first investigate the accuracy achieved with the original method as a function of the number of nodes N and the non-linear HOS order M . We choose a case where ka = 0.1 and kh = 0.92. This case presents a small relative water depth and a medium steepness. Indeed, the highest stable wave is given by [START_REF] Williams | Limiting gravity waves in water of finite depth[END_REF] and the parametrized formula of [START_REF] Fenton | Nonlinear wave theories[END_REF]. In this case of shallow water conditions (kh = 0.92), the value of the steepness (ka = 0.1) represents 33% of the maximum curve.

Fig. 5 presents the logarithm of the relative error made on the vertical velocity as a function of N and M for two values of β h0 . The lower surface is for β h0 = 0% (i.e the initial HOS method) while the upper surface represents the convergence for a variation of the bottom β h0 = 25%. We observe that when β h0 increases, the error on the vertical velocity becomes higher than the error with a flat bottom. Nevertheless, the method still converges exponentially and preserves the accuracy for perturbations of the bottom equal to 25%. In this particular case, one can notice that the convergence is more difficult to achieve in the Mdirection, because of the large non-linearities induced by the bottom. The order of non-linearity at the bottom probably needs to be larger than the one at the free-surface, showing a possible advantage of the improved method presented hereafter.

We can compute cases where β h0 = 50%. The results are presented in Fig. 6. We still observe the exponential convergence of the method, but with a reduced rate compared to the previous case (as we saw just before, the convergence is We conclude that with the original method, when kh increases (with ka kh=0.916

β h0
k.a=0.01 k.a=0.05 k.a=0.1 0% 1.7 10 -14 4.0 10 -10 2.5 10 -7 5% 4.8 10 -12 1.8 10 -9 3.0 10 -7 25% 5.5 10 -8

2.6 10 -6 3.4 10 -5 50% 3.9 10 -6

1.1 10 -4 1.2 10 -3 ka=0.1

β h0
k.h=0.67 k.h=3 k.h=10 0% 5.7 10 -8 6.3 10 -9 6.2 10 -9 5% 2.2 10 -7 6.3 10 -9 6.2 10 -9 25% 3.9 10 -5

1.8 10 -6 2.6 10 -7 50% 1.1 10 -3

3.8 10 -4 1.8 10 -4

Table 1: Relative error on W , N = 64, M = 8.

and β h0 fixed) the convergence rate is better. This result was expected because when h λ > 2 (i.e. kh > π), we may consider to be in a deep water problem with no influence of the bottom on wave propagation. Thus the bottom has less influence and the convergence rate is higher. When ka and β h0 increase (with the other parameters fixed), the convergence rate decreases. Of course there is a limit in the admissible values of ka and β h0 . The steepness can not exceed the theoretical highest stable wave given by the formula of [START_REF] Fenton | Nonlinear wave theories[END_REF] and β h0 has to remain lower than 100%. So this test case (which is very demanding) shows that a quite large range of parameters can be treated with the original HOS method.

Improved method

The same test case is studied with the improved method. Remember that the influence of the order of non-linearity of the bottom M b has to be investigated in addition to other parameters. We thus fix N = 32 (we have checked that the method converges with N with an exponential convergence rate), and we observe the influence of M and M b for different values of the parameters kh, ka, and β h0 . We observed with the original method that the convergence was more difficult to achieve in the M -direction for the case ka = 0.1, kh = 0.92 and β h0 = 25%. We guess that with the improved method we have to increase M b (more than M ) in order to reach the convergence. The error achieved on the vertical velocity is presented in Fig. 7. In this figure, we can clearly see that the method converges exponentially with M and M b , which is a very interesting feature of the new improved scheme. Of course the method also converges exponentially with N but it can not be seen on the figure. The same conclusions as before can been drawn with increasing kh, ka and β h0 and comparisons between the two approaches are provided in the following section.

Comparison of the two methods

As both methods present an exponential convergence rate (with N and with the different orders of non-linearity M and M b ), it is now interesting to investigate which method is more accurate, depending on the values of kh, ka, and β h0 . (ka = 0.01) shows that at a fixed value M , we can always find a M b (threshold) to get at least the same error on the vertical velocity as with the original method.

Moreover, if M b is higher than this threshold value, the accuracy is better. Thus the improved method can improve the accuracy achieved on the vertical velocity.

For steeper cases, Fig. 9, (ka = 0.15) proves that the improved method is always more accurate than the original one (if M b > 1). It means, in the case of Fig. 9, that the non-linearities due to the free-surface are more important than in the case of Fig. 8, because of the higher steepness. The influence of β h0 is represented on Fig. 10 and Fig. 11. In Fig. 10, for M b ≥ 3 the error achieved on the vertical velocity is converged so the corresponding curves are superimposed. For theses cases we see that for a small variation of the bottom (Fig. 10), both methods present the same accuracy. But when β h0 increases, the improved method seems much more accurate, for the same reason as before: when β h0 is higher, high non-linearities induced by the bottom are more important than those due to the free-surface, so a larger M b is needed. Fig. 11 also shows the capacity of the improved scheme to treat high variations of the bottom (here β h0 = 75% with a medium steepness ka = 0.15 and a small relative water depth kh = 0.916).

The influence of the relative water depth kh can be seen on Fig. 12 and Fig. 13.

In Fig. 12 the curves corresponding to M b ≥ 3 are superimposed. For small kh both methods show almost the same accuracy because the non-linearities on the free-surface elevation are important, and thus the accuracy is governed by the order M . However, for higher kh, the improved method allows a better accuracy, because as before, the non-linearities induced by the bottom become more important than those due to the free-surface.

An important conclusion from these investigations is that, with the improved method, we can always find numerical parameters to achieve at least the same accuracy as with the original method. Thus the improved method is more accurate if the different orders of non-linearity are correctly chosen. Moreover, in some cases, the use of a couple (M, M b ) rather than a couple (M, M ) should reduce the computational effort, while maintaining the accuracy (see 5.1 below).

Once the vertical velocity is correctly computed, one can advance the solution in time. One way to estimate whether the propagation is correct is to check the phase shift with respect to the reference numerical solution of [START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF]. A phase shift is usually observed after a long propagation time with any numerical method. The HOS method is very accurate with respect to this phase shift (see [START_REF] Bonnefoy | Time-domain simulation of nonlinear water waves using spectral methods[END_REF]). It happens even if the bottom is flat (ie β h0 = 0%), and we have checked in [START_REF] Gouin | Development and validation of a highly nonlinear model for wave propagation over a variable bathymetry[END_REF] that neither method induces a prohibitive phase shift.

Computational effort

The computational effort is directly proportional to the number of FFTs used in the HOS procedure. For the initial HOS method with a flat bottom (Eqs. 10 and 11), the number of FFTs required is M M +3 2 . To account for a variable bathymetry, the additional number of FFTs used in the original method (Eq.16) is:

n = 2M (M -1) (26)
While for the improved method this number is:

n = 4M b M (27)
Even if the additional cost induced by a variable bathymetry is in O(M 2 ), we are confident in the efficiency of the two HOS methods which are very interesting in terms of CPU time compared to other existing non-linear potential flow solvers.

Moreover, we observe that if M = M b , the original method will be more efficient. On the contrary, the improved method should be more adequate if:

M b < M -1 2 or M < M b -1 2 (28)
As we saw just before, the improved method allows to achieve errors on the vertical velocity less than or equal to those obtained with the original method, depending on the value of M b . Thus, in some cases where the two needed orders of non-linearities are very different, the use of a couple (M, M b ) rather than a couple (M, M ) will reduce the computational effort (while keeping the accuracy).

In order to confirm the asymptotic scaling of the computational effort

N d log(N d )
of the proposed models, Fig. 14 presents the computational time for 100 periods with respect to the number of dealiased points N d in a computational domain of 50 wavelengths. Here

N d = N M +1 2
as we perform a total dealiasing. This figure was obtained with a HOS order fixed to M = 8 and a variation of the bottom β h0 = 40%. The initial HOS method with a flat bottom and M = 8 is also presented for comparison. The dashed lines represent the expected linear scaling in N d log(N d ). We observe that the proposed models and the initial HOS method scale the same way with N d log(N d ), which is a very interesting feature of the new HOS methods. Moreover, with this choice of parameters (M = 8 and β h0 = 40%), the computational time needed to compute a variable bathymetry is three times the computational effort necessary for a flat bottom. Looking at the number of FFTs required, the initial method for M = 8 needs n = 44 FFTs while the original method requires n = 122 + 44 = 156 FFTs (cf Eq.26). The theoretical ratio is thus 3.5 which agrees well with the empirical result of 3.

Partial conclusion

Both methods can compute various cases depending on the value of β h0 , kh and ka. They present a good accuracy as well as an attractive scaling of the computational effort with increasing problem size. Moreover, for increasing steepness, the error on the vertical velocity becomes larger and the convergence is slower, but still exists. The improved method is more flexible as we can dissociate the different orders of non-linearity. Furthermore, we can always find numerical parameters to achieve a better accuracy than with the original method. The computational effort, which scales as N d log(N d ) (with N d the number of dealiased points), will also probably be reduced in some configurations if the orders of non-linearity needed on the free-surface and on the bottom are different. We will compare the required computational time in the applications cases given in the following.

Thus this highly demanding test case demonstrates the ability of our methods to represent non-negligible variations of the bottom with high accuracy and efficiency. We are now confident on the applicability of these schemes to realistic problems.

Validation

Two test cases have been treated to demonstrate the ability of the proposed schemes to model wave propagation over a variable bathymetry in realistic configurations. As an example of a small bottom variation around a mean water depth, the first test case considers Bragg reflection from a sinusoidal bottom patch. The second one intends to focus on a larger bottom variation with the shoaling phenomenon.

Bragg reflection

This case aims at representing a small bottom variation around a mean water depth in order to satisfy the conditions of the Taylor expansion (see Eq.13). We consider the development of incident and Bragg reflected waves over a sinusoidal bottom patch of wavenumber k b (see Fig. 15). If the class I Bragg condition is satisfied, the reflected wave should be amplified as a result of resonant quadratic interaction between the incident wave and the bottom variation. This case has been studied experimentally by [START_REF] Davies | Surface-wave propagation over sinusoidally varying topography[END_REF] and numerically by [START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF]; [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF]; [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF].

For small incident waves and small bottom slopes, reflection near Bragg resonances is well predicted by the multiple-scale perturbation theory of [START_REF] Mei | Resonant reflection of surface water waves by periodic sandbars[END_REF]. Here we analyse non-linear effects using the original method. Note that the improved method gives exactly the same results. The conditions of [START_REF] Davies | Surface-wave propagation over sinusoidally varying topography[END_REF] are used to compare with experiments. The bottom ripple patch is defined as the variation around the mean water depth h = h 0 + β(x) by :

β(x) = d sin(k b x) for -L 0 2 ≤ x ≤ L 0 2 (29)
as depicted in Fig. 15. The free surface is located at z = 0m and the mean water depth is at h 0 = 1m. The ripple slope is k b d = 0.31, the ripple amplitude is d = 0.1m and the length of the patch is L0 λ b = 10 (i.e. we consider a patch of 10 sinusoidal ripples of wavelength λ b = 2π k b ). The incident wave is set at the linear resonance condition of k = k b 2 with a wave steepness ka = 0.05. We perform simulations with N λ = 16 nodes per wavelength and a HOS order M = 2 with the original method (the equivalence is obtained with M = 2 and M b = 1 for the improved method). The steady-state is reached after t = 41T . This choice of order of non-linearity on the free-surface and on the bottom variation is sufficient to obtain converged results, since main physical phenomena are accounted for by Class I Bragg reflection which is second order, as explained in [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF]. The local reflection coefficient R(x) is then extracted using the method of [START_REF] Suh | Separation of incident and reflected waves in wavecurrent flumes[END_REF] and analysed hereafter.

The results are presented in Fig. 16 along with the experimental measurements of [START_REF] Davies | Surface-wave propagation over sinusoidally varying topography[END_REF] and the solution given by the perturbation theory of [START_REF] Mei | Resonant reflection of surface water waves by periodic sandbars[END_REF]. Ripples are located in x/λ b ∈ [-5; 5]. It appears clearly that our numerical results accurately reproduce the experiments and are in excellent agreement with the perturbation theory. Both methods (original and improved) give exactly the same results as those presented in [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF]. Thus the proposed methods are validated for a classical configuration with a small variation of the bottom ( β h0 = 0.1). We will now focus on the next test cases with higher bottom variations, to check the ability of the methods to treat realistic bathymetry profiles with non-negligible variations.

Shoaling of linear waves

Here we consider the shoaling of linear waves which are deformed during their propagation from deep to shallow water (see Fig. 17). This validation case has been run with the improved method only. Indeed, it requires a linear solution at the free-surface (to compare our results to an analytical solution) and a highly non-linear order at the bottom to account for the influence of the Taylor expansion around a mean water depth. This test case has been studied for instance in [START_REF] Bingham | A Fourier-Boussinesq method for nonlinear water waves[END_REF]; [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF].

We used a bottom profile defined by: with h 0 = 0.5 m, h 1 = 0.1 m and L = 24 m. Here β h0 = 0.8 so the hypothesis of the Taylor expansion Eq.17 is not fulfilled. Nevertheless, we want to investigate if the method can deal with a high relative variation of the bathymetry. The bottom topography and the incident conditions allow a transition from deep to shallow water (from kh = 1.57 to kh = 0.56). A linear wave of amplitude a 0 = 5.10 -3 m and wavelength λ 0 = 2m is generated at the left end of the domain and is propagated until a steady-state is reached everywhere. A relaxation zone at the left of the domain generates the incoming wave and another zone at the right of the domain absorbs the waves after shoaling to prevent reflection and to ensure periodicity. Relaxation zones are thus composed of one absorbing zone (5 wavelengths), one generation zone (3 wavelengths), the observation zone (12 wavelengths) and a second absorbing zone (10 wavelengths) (see Fig. 17).

β(x) = h 0 -h 1 2 1 + tanh sin (πx/L) 1 -(2x/L) 2 , - L 2 ≤ x ≤ L 2 (30)
The mean water depth h 0 is chosen at 0.3 m to reduce the maximum β h0 (from 0.8 to 0.66 at its maximum value) and to ensure a convergence with respect to M b .

The steady-state is reached after t = 27T . Fig. 18 shows the computed steady state amplitude envelope (M = 1 and M b = 15 with h 0 = 0.3m) and the one predicted by energy conservation from linear theory: We see that the amplitude first decreases slightly with the bottom slope, and then increases and stabilizes. We note that our model compares accurately with the linear prediction. The increasing value of the wavenumber as the wave travels from deep to shallow water is also an expected phenomena (see [START_REF] Bingham | A Fourier-Boussinesq method for nonlinear water waves[END_REF]). Thus this test case validates the improved method even for non-negligible relative bottom variations β h0 .

a a 0 = k (1 + 2h 0 k 0 / sinh(2h 0 k 0 )) k 0 (1 + 2hk/ sinh(2hk)) 1/2 (

Application : Harmonic Generation over a Submerged Bar

These application cases consider the transformation of a non-linear monochromatic wave as it travels up and over a submerged bar. As the waves travel up the front slope of the bar, they steepen dramatically, and higher harmonics are generated due to non-linear interactions. These higher harmonics produce an irregular pattern behind the bar. These validation cases have often been used as a discriminating test case for non-linear models of surface wave propagation over a variable bottom because they require an accurate propagation of waves in both deep and shallow water. The first case presented considers a bar with mild slopes while the second one focuses on steeper slopes.

Mild slope

This test case has been studied both experimentally by [START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF]; [START_REF] Luth | Kinematics of waves breaking partially on an offshore bar[END_REF] and numerically by [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF]; [START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF]. The bottom variation is defined by:

β (x) =                0.
and can be seen in Fig. 19. It has been scaled with a factor of two in comparison with the experiment of [START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF]. The value of the front and back slopes are respectively 5% and 10%, corresponding to classical orders of magnitude observed on continental shelves. The relative variation of the bottom is β h0 = 0.75 and we want to investigate if the method can deal with such variations.

Regular waves are generated at the left side of the domain using the solution of [START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF] of steepness ka = 0.017 and relative water depth kh = 0.67. The period is fixed to 2.02s with an amplitude of 0.01m. With the two methods the convergence and steady-state are reached after t = 25T with a number of nodes per wavelength N λ = 40. The original method needs an HOS order M = 17 to reach the steady-state. Indeed, such a high-order is needed on the bottom because non-linearities on the free-surface are small (see the value of ka) while the ones induced by the bottom variation are important.

As a comparison, the improved method only needs an order M = 3 on the free surface and an order M b = 16 on the bottom to reach the steady-state.

When the convergence and steady-state are achieved, exactly the same results are obtained with the two methods. Thus, in the following, we present only the results with the improved method. Moreover, we demonstrated in that when the two needed orders of non-linearities are very different, the use of a couple (M, M b ) rather than a couple (M, M ) will reduce the computational effort. Thus, in this case, the use of the improved method rather than the original method improves the efficiency, while keeping the accuracy.

Time series of free surface elevations

For the incident wave conditions (T 0 , a 0 )= (2.02s, 0.01m), a snapshot of the surface elevation (scaled by a factor of 3) is presented in Fig. 20 and The comparison between our numerical results and the experimental data is very good, and similar to the results obtained with other numerical methods [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF]; [START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF]). In particular, the asymmetry of the shoaling waves and the generation of higher harmonics are correctly reproduced by the two methods. Thus, both free-surface non-linearities and bottom non-linearities are correctly solved and we are confident in the accuracy of the model and its ability to treat large bottom variations.

Harmonic analysis

For a deeper comparison, a harmonic analysis of the surface elevation is run as presented in [START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF]. Our results are shown in Fig. 22 where the normalized amplitude of the first five harmonics of the wavefield are depicted as functions of x. As expected, we can clearly observe the generation of high-harmonics over the bar. Moreover, the comparison of all harmonics with the measurements is good even up to the fifth-harmonic, indicating that non-linear effects due to the bottom variation are correctly resolved. It is also very similar to the numerical results presented in [START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF], especially with similar slight discrepancies observed on the first harmonic. The oscillations are due to the interaction between the incident waves and the waves reflected on the first slope of the bar. We can notice that even if the submerged bar (which represents a very large relative bottom variation), the steepness is weak (ka = 0.017), and the bottom variation does not take place in the whole domain, helping the model to solve this problem accurately.

Looking at the computational effort, we observe that the improved method is more efficient: the CPU time needed to compute one period is 5.5s with the improved method (M = 3 and M b = 16) and 8.6s with the original method (M = 17).

Steep slope

We now apply the two HOS methods to a submerged bar with steeper front and back slopes, as presented in the experimental set-up of [START_REF] Ohyama | Applicability of numerical models to nonlinear dispersive waves[END_REF].

In this case, β h0 = 0.7 and the value of the slope is 50% so it allows to check the domain of validity of our model. This experiment has also been used by [START_REF] Shen | Simulation of wave propagation over a submerged bar using the VOF method with a two-equation k-turbulence modeling[END_REF]; [START_REF] Gobbi | Wave evolution over submerged sills: Tests of a high-order Boussinesq model[END_REF]; [START_REF] Benxia | Wave decomposition phenomenon and spectrum evolution over submerged bars[END_REF] as a validation case for their numerical methods. of the surface elevations at locations a (over the bar) and b (after the bar) (see Fig. 23) are shown in Fig. 24. The experimental data come from the experiment of [START_REF] Ohyama | Applicability of numerical models to nonlinear dispersive waves[END_REF]. The comparison between our numerical results and the experimental data is also good even if we observe slight discrepancies with respect to the experiments. In this steeper case, the maximum variation of the bottom is almost the same as in the mild slope case ( β h0 = 70%) but the gradient is much larger. Nevertheless, the large slopes are more difficult to simulate with the proposed 570 models. However, in this extreme configuration, we still manage to obtain reasonable accuracy on the free-surface elevation. Moreover, the level of accuracy achieved on the time series of the free-surface elevations is better than the one (where the two orders of non-linearities are very different), the use of a couple (M, M b ) rather than a couple (M, M ) will even reduce the computational effort (while keeping or improving the accuracy). The flexibility and efficiency of the improved method are thus important features of the new scheme.

Then, by a series of two validation cases, we have shown the accuracy and efficiency of the two methods. Efficient procedures for waves generation/absorption have been developed to ensure the periodicity of the simulation. The first test case reproduces Bragg reflection over small bottom variations and shows results in agreement with [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF]; [START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF]; [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF]. The second test case was implemented to reproduce the shoaling of linear waves. It shows very good agreement with the analytical solution and with numerical results described in [START_REF] Bingham | A Fourier-Boussinesq method for nonlinear water waves[END_REF]; [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF].

Finally, the application case of a submerged bar which simulates large and realistic bottom variations has been presented. It shows very good agreement with both numerical and experimental data, and thus proves the ability of the two methods to accurately simulate the propagation of waves over large variations of the bathymetry and of its gradient. The improved method seems to be more efficient on the specific case of the underwater bar because the bottom and the free-surface do not require expansions with the same order of non-linearity, so the computational time is reduced with this method. Thus, we have demonstrated the ability of our methods to compute high water depth variations and bottom slopes. Despite the limitations inherent to the Taylor expansions, very good results have been obtained even for very steep waves over highly varying bathymetries. Moreover, the two iterative methods presented appear very effective, and the choice of one or another will depend on the computed case. They both can be used to accurately characterize the wavefields in coastal regions. The calculation of corresponding wave loads on marine renewable energy systems is the next step. These loads may be directly computed from wave kinematics using added mass and drag terms estimated His PhD was dealing with the modeling of non-linear processes during generation and propagation of gravity waves.

Figure 1 :

 1 Figure 1: Description of the fluid domain

Figure 2 :

 2 Figure 2: Use of the bottom condition in the temporal solution.

Figure 3 :

 3 Figure 3: Temporal solution of the FSBC.

  from Liu and Yue (1998)) is depicted in Fig.2and explained hereafter.

  l = 0 to M b (with M b the order of non-linearity of the bottom):

Figure 4 :

 4 Figure 4: Schematic view of the validation case

Figure 5 :Figure 6 :

 56 Figure 5: Error on the vertical velocity in function of N and M -kh=0.916 and ka=0.1. Lower-surface : β h 0 = 0%. Higher-surface : β h 0 = 25%.

Figure 7 :

 7 Figure 7: Error on the vertical velocity in function of M and M b -kh = 0.916, ka = 0.1 and β h 0 = 25%.

Fig. 8

 8 Fig.8 and Fig.9 present the error achieved on the vertical velocity as function of M and M b for β h0 = 40%, kh = 0.92 and two different values of ka. In Fig.9, from M b = 10 the error achieved on the vertical velocity is converged so the three corresponding curves (M b = 10, M b = 15 and M b = 20) are superimposed. Fig.8

Fig. 8

 8 Fig.8 also highlights another distinction between the two methods. Indeed, one could think that the original method with M = 10 and the improved method with M = 10 and M b = 10 would lead to the same results. Nevertheless, when computing the modal amplitudes B

  at each step m, the improved method calculates more terms (if M b > M ), inducing a difference in the error achieved on the vertical velocity.

Figure 8 :

 8 Figure 8: Error on the vertical velocity as function of M and M b -β h 0 = 40% -kh = 0.916 and ka = 0.01.

Figure 9 :

 9 Figure 9: Error on the vertical velocity as function of M and M b -β h 0 = 40% -kh = 0.916 and ka = 0.15.

Figure 10 :

 10 Figure 10: Error on the vertical velocity as function ofM and M b -β h 0 = 10% -kh = 0.916 and ka = 0.15.

Figure 11 :

 11 Figure 11: Error on the vertical velocity as function of M and M b -β h 0 = 75% -kh = 0.916 and ka = 0.15.

Figure 12 :

 12 Figure 12: Error on the vertical velocity as function of M and M b -β h 0 = 40% -kh = 0.5 and ka = 0.1.

Figure 13 :

 13 Figure 13: Error on the vertical velocity as function of M and M b -β h 0 = 40% -kh = 3 and ka = 0.1.

Figure 14 :

 14 Figure 14: Computational effort with respect to N . M = 8.

Figure 15 :

 15 Figure 15: Bottom topography with a patch of 10 sinusoidal ripples of amplitude d = 0.1 and slope k b d = 0.31.

Figure 16 :

 16 Figure 16: Bragg reflection from a sinusoidal bottom ripple patch over -5λ b ≤ x ≤ 5λ b with ka = 0.05 and k b d = 0.31.

Figure 17 :

 17 Figure 17: Computed wave profile (scaled by a factor of 10) of a linear wave of amplitude a 0 = 5.10 -3 m and wavelength λ = 2m over a shoaling bottom profile at steady-state.

Figure 18 :

 18 Figure 18: Comparison between the envelope predicted by the linear theory and the envelope obtained with the computed wave profile.

Figure 19 :

 19 Figure19: Sketch of the submerged bar as described in the experiments of[START_REF] Dingemans | Comparison of computations with Boussinesq-like models and laboratory measurements[END_REF]. Blue markers correspond to the gauges in the experiments.

  the time histories of the surface elevations at various locations are shown in Fig.21. The location of the gauges can be seen in Fig.19. The experimental data comes from the experiment of Luth et al. (1994).

Figure 20 :

 20 Figure 20: Snapshot of the surface elevation at a steady-state (scaled by a factor of 3).

Figure 21 :

 21 Figure 21: Time series of measured and computed surface elevations at 6 positions. Solid line: calculations. Points: measurements Luth et al. (1994).

Figure 22 :

 22 Figure 22: Harmonic analysis at steady-state for the mild slope case. N λ = 40, M = 3 and M b = 16.

Figure 24 :

 24 Figure 24: Time series of measured and computed surface elevations at 2 positions. Solid line: calculations. Points: measurements Ohyama et al. (1995).

765

  Pierre Ferrant is professor of hydrodynamics in ECN, and is managing the LHEEA Lab. of ECN. His research covers free surface flows, ocean waves modelling and wave structure interactions, using both numerical and experimental approaches. He graduated in Ecole Centrale de Nantes in 1981, obtained a Master in Ocean engineering in 1982, and a PhD in 1988. After a 10-year period in 770 a private research company, he began an academic career in ECN in 1998. Over the last 15 years, Pierre Ferrant has supervised 25 PhD thesis and published some 200 papers in refereed journals and international conferences.

https://github.com/LHEEA/HOS-ocean/wiki

The bottom variation is presented in Fig. 23 and defined for this steeper case by: The initialisation of η and φ is based on a solution of [START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF] of steepness ka = 0.039 and relative water depth kh = 0.77. The period is fixed to 2.01s with an amplitude of 0.025m. The relaxation zones are the same as previously.

Simulations are run until the steady-state is reached after t = 30T with a number of nodes per wavelength N λ = 60 and an HOS order M = 17. The improved method only needs an order M = 5 on the free surface and an order M b = 16 on the bottom to reach the steady-state. As previously, the results are identical with the two methods, so we only present the results achieved with the improved one. The computational effort of the two methods is different because it depends from the numerical parameters M and M b , as seen before.

For the incident wave conditions (T 0 , a 0 )= (2.01s, 0.025m), the time histories presented in [START_REF] Shen | Simulation of wave propagation over a submerged bar using the VOF method with a two-equation k-turbulence modeling[END_REF] and almost the same as in [START_REF] Gobbi | Wave evolution over submerged sills: Tests of a high-order Boussinesq model[END_REF]; [START_REF] Benxia | Wave decomposition phenomenon and spectrum evolution over submerged bars[END_REF].

With the two proposed methods, the bottom and the gradient need to be continuous, and we clearly see that this is not the case in theses two configurations. Nevertheless, we reach to treat the two cases of the underwater bars with mild and steep slopes. For more complicated bathymetries, the bathymetry could be smoothed to avoid the development of instabilities.

Thus, these HOS models appear suitable to simulate wave propagation with large variations of water depth and relatively large bottom slopes.

Conclusion

In the frame of the initial HOS model developed for a flat bottom, we have implemented two numerical methods to simulate non-linear free surface waves over variable depth. Both methods are based on a Taylor expansion of the bottom boundary condition with respect to the mean water depth. The first method uses the same order of non-linearity for the bottom and for the free-surface as presented in [START_REF] Liu | On generalized Bragg scattering of surface waves by bottom ripples[END_REF], while the second method is an extension of the work of [START_REF] Guyenne | A high-order spectral method for nonlinear water waves over moving bottom topography[END_REF] to the HOS formalism by considering two independent orders of non-linearity.

We first characterized and demonstrated the validity of both methods by computing large but constant bottom variations. We demonstrated that the error made on the vertical velocity converged exponentially with increasing N , M and M b . The improved method is more flexible as we can dissociate the different orders of non-linearity. Moreover, the computational effort to account for a variable bathymetry scales with N d log(N d ) for both methods (where N d is the number of dealiased nodes), similarly to the initial HOS method. We also saw that we can always find numerical parameters (M, M b ) leading to a better accuracy with the improved method than with the original one. In some cases from Morison's formula. A more reliable and accurate solution will be to couple one of the proposed models for wave propagation, with a CFD solver accounting for the interaction with the structure.

The extension of the model to irregular waves and to 3D cases is conceptually straightforward with the use of relaxation zones and is under development.