Cécile Raoult
email: cecile.raoult@gmail.com

Michel Benoit
email: benoit@irphe.univ-mrs.fr

Marissa L Yates
email: marissa.yates-michelin@cerema.fr

Development and validation of a 3D RBF-spectral model for coastal wave simulation

Keywords: nonlinear, dispersive, water waves, potential theory, Zakharov equations, Radial Basis Functions

 is extended to three dimensions (3D). To maintain the flexibility of the approach with the goal of applying the model to irregularly shaped domains, the horizontal plane is discretized with scattered nodes. The horizontal derivatives are then estimated using the Radial Basis Function -Finite Difference (RBF-FD) method, while a spectral approach is used in the vertical dimension. A sensitivity analysis examined the robustness of the RBF-FD approach as a function of RBF parameters when estimating the derivatives of a representative function. For a targeted stencil size between 20 and 30 nodes, Piecewise-Smooth (PS) polyharmonic spline (PHS) functions are recommended, avoiding the use of Infinitely-Smooth (IS) RBFs, which are less appropriate for the desired applications because of their dependence on a shape parameter. Comparisons of simulation results to observations from two wave basin experiments show that nonlinear effects induced by complex bottom bathymetries

Introduction

Accurate wave propagation models are required for a wide range of coastal management and engineering applications, including the design of coastal structures and the evaluation of coastal risks. In the nearshore zone, nonlinear and dispersive effects, characterized by the wave steepness and relative water depth, respectively, can be particularly important. The wave steepness = kH/2 and relative wave height H/h (where H, k, and h are the characteristic local wave height, wave number, and water depth, respectively) increase significantly as waves shoal in shallow water and approach the breaking point. The relative water depth µ = kh is often large in deep and intermediate water, or for the shorter waves in the sea state, indicating the importance of dispersive effects. Wave propagation models thus need to capture properly nonlinear and dispersive effects to simulate accurately offshore and coastal engineering problems. Two-dimensional cross-shore (2DV) wave models can be used as a preliminary step in coastal studies, but 3D models are needed to capture fully the effects of alongshore bathymetric variations, variable wave incidence, the presence of coastal or harbor structures, etc.

A wide variety of mathematical models exist to simulate nearshore wave propagation, and [START_REF] Yates | Modélisation non-linéaire et dispersive des vagues en zone côtière: étude comparative de deux méthodes de simulation précises[END_REF] and [START_REF] Benoit | A comparison of simulation approaches based on the Zakharov equations for nonlinear waves in the coastal zone[END_REF] chose to develop a 2DV fully nonlinear and dispersive potential flow model based on the Euler-Zakharov equations as a compromise between accuracy, mathematical complexity, and computational time. A variety of other approaches exist, and recent non-exhaustive reviews are summarized by [START_REF] Raoult | Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments[END_REF] and [START_REF] Gouin | Development and validation of a non-linear spectral model for water waves over variable depth[END_REF], for example. The chosen approach, based on fully nonlinear potential flow (FNPF) theory, requires solving the Laplace boundary value problem (BVP), which is implemented numerically in the model using a spectral approach in the vertical direction [START_REF] Tian | A numerical model on the interaction between nearshore nonlinear waves and strong currents[END_REF] and finite difference schemes in the horizontal direction. The nonlinear and dispersive capacities of the 2DV version of the model were demonstrated in Yates andBenoit (2012, 2015) and [START_REF] Raoult | Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments[END_REF]. The accuracy of a similar approach using the Chebyshev-Tau method in the vertical and a Fourier collocation method in the horizontal was demonstrated in [START_REF] Christiansen | Efficient hybrid-spectral model for fully nonlinear numerical wave tank[END_REF].

Extending 2DV modeling approaches to 3D increases significantly the number of numerical challenges to overcome, including but not limited to the computational time, domain geometry, and boundary condition specification or far-field representation. Therefore, a variety of different numerical approaches have been used to solve the 3D FNPF problem by reducing the dimension of the problem (e.g. boundary element [START_REF] Romate | Boundary integral equation formulations for free-surface flow problems in two and three dimensions[END_REF][START_REF] Grilli | A fully non-linear model for threedimensional overturning waves over an arbitrary bottom[END_REF][START_REF] Fochesato | Numerical modeling of extreme rogue waves generated by directional energy focusing[END_REF][START_REF] Nimmala | An Efficient Three-Dimensional FNPF Numerical Wave Tank for Large-Scale Wave Basin Experiment Simulation[END_REF] or High-Order Spectral (HOS) methods [START_REF] Ducrozet | 3-D HOS simulations of extreme waves in open seas[END_REF][START_REF] Ducrozet | A modified High-Order Spectral method for wavemaker modeling in a numerical wave tank[END_REF]), discretizing the whole domain but taking local derivatives only (e.g. finite element [START_REF] Wu | Numerical simulation of sloshing waves in a 3D tank based on a finite element method[END_REF] or finite difference schemes [START_REF] Engsig-Karup | An efficient flexibleorder model for 3D nonlinear water waves[END_REF]), or searching for fast numerical methods for solving the Laplace problem (e.g. integral equations [START_REF] Fructus | An efficient model for three-dimensional surface wave simulations: Part I: Free space problems[END_REF] or coupled modes [START_REF] Belibassakis | Directional wave spectrum transformation in the presence of strong depth and current inhomogeneities by means of coupled-mode model[END_REF]), each approach having its own advantages and disadvantages.

A number of 3D Numerical Wave Tanks (NWTs) have been developed using high-order Boundary Element Methods (BEM) [START_REF] Romate | Boundary integral equation formulations for free-surface flow problems in two and three dimensions[END_REF][START_REF] Grilli | A fully non-linear model for threedimensional overturning waves over an arbitrary bottom[END_REF][START_REF] Fochesato | Numerical modeling of extreme rogue waves generated by directional energy focusing[END_REF][START_REF] Nimmala | An Efficient Three-Dimensional FNPF Numerical Wave Tank for Large-Scale Wave Basin Experiment Simulation[END_REF], which are an efficient and accurate approach for solving the Laplace BVP by reducing the dimensionality of the discretized problem. However, standard BEM techniques yield nonsymmetric and fully populated matrices that require fast solution methods (e.g. fast multipole algorithm [START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF]) or advanced numerical implementations (e.g. pre-corrected Fast Fourier Transform methods [START_REF] Yan | An efficient high-order boundary element method for nonlinear wavewave and wave-body interactions[END_REF]) to avoid becoming computationally prohibitively expensive to solve. This is especially the case for simulating irregular nonlinear waves that span a wide range of wavelengths, thus requiring fine grids.

In the literature, two other methods have been used to extend 2DV approaches to 3D with Finite Element Methods (FEM) (e.g. 2DV [START_REF] Wu | Finite element analysis of two-dimensional non-linear transient water waves[END_REF] to 3D [START_REF] Wu | Numerical simulation of sloshing waves in a 3D tank based on a finite element method[END_REF]) or finite difference schemes (e.g. 2DV [START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF] to 3D [START_REF] Engsig-Karup | An efficient flexibleorder model for 3D nonlinear water waves[END_REF]). These approaches require discretizing the entire fluid domain, but the only non-zero elements in the coefficient matrix are the neighboring points (with the number of points depending on the chosen order of the numerical scheme). While the full 3D Laplace problem still must be solved at each time step, recent work has improved the numerical efficiency of such codes (e.g. using preconditioned defect correction methods in 2DV [START_REF] Christiansen | Efficient hybrid-spectral model for fully nonlinear numerical wave tank[END_REF] or in 3D in the OceanWave3D code (Engsig-Karup, 2014)) that may then be used for coastal and offshore engineering applications.

Another even more computationally efficient approach uses the HOS method (e.g. [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF]; [START_REF] West | A new numerical method for surface hydrodynamics[END_REF]; [START_REF] Craig | Numerical simulation of gravity waves[END_REF]), which is based on a Taylor series expansion of the velocity potential about a reference water level (often the mean water line). This approach was originally developed for unbounded domains and then for fi-nite depth cases, with periodic lateral boundary conditions. [START_REF] Ducrozet | 3-D HOS simulations of extreme waves in open seas[END_REF][START_REF] Ducrozet | A modified High-Order Spectral method for wavemaker modeling in a numerical wave tank[END_REF] developed an open-source model, HOS-ocean, which is an extension of this method to take into account the generation and propagation of regular and irregular, multidirectional waves. Recently, HOS-ocean was also extended to simulate wave propagation over variable bottom bathymetries by [START_REF] Gouin | Development and validation of a non-linear spectral model for water waves over variable depth[END_REF].

Additional approaches for solving the BVP include a fast Laplace equation solver using integral equations and an iterative solution procedure that converges rapidly (in 2DV [START_REF] Clamond | A fast method for fully nonlinear water-wave computations[END_REF] and then 3D [START_REF] Fructus | An efficient model for three-dimensional surface wave simulations: Part I: Free space problems[END_REF]), or a fully dispersive coupled-mode model (in 2DV [START_REF] Belibassakis | A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions[END_REF] and then 3D [START_REF] Belibassakis | Directional wave spectrum transformation in the presence of strong depth and current inhomogeneities by means of coupled-mode model[END_REF]).

Here, the objective is to extend the two dimensional model developed by [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF] and [START_REF] Raoult | Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments[END_REF] to three dimensions, maintaining the efficient and accurate spectral approach used in the vertical. The 2DV code used finite difference schemes to estimate horizontal derivatives. This approach is still applicable in two horizontal dimensions for simple domain geometries that can be discretized by regular meshes. However, the model cannot then be applied easily to complex domains. Therefore, to overcome these limitations, a meshless approach based on Radial Basis Functions (RBF) is implemented in the code and tested to propagate highly nonlinear and dispersive waves. The RBF method has been used extensively in other fields of research for a variety of physical and engineering problems, including, for example, diffusion [START_REF] Sarler | Meshfree explicit local radial basis function collocation method for diffusion problems[END_REF][START_REF] Chen | Fractional diffusion equations by the kansa method[END_REF], radiative transport [START_REF] Kindelan | Application of the RBF meshless method to the solution of the radiative transport equation[END_REF], combustion [START_REF] Kansa | Numerical simulation of twodimensional combustion using mesh-free methods[END_REF], shallow-water models [START_REF] Zhou | A grid-free, nonlinear shallow-water model with moving boundary[END_REF][START_REF] Hon | Multiquadric solution for shallow water equations[END_REF], and flow simulations using the Navier-Stokes equations [START_REF] Shu | Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier-Stokes equations[END_REF]. Here, the RBF method is applied to evaluate its capacity to simulate accurately nonlinear wave propagation.

This paper presents a brief literature review of RBF methods (Section 2), before describing the mathematical model and its numerical implementation in the code (Section 3), focusing on the RBF Finite Difference (RBF-FD) method used to estimate the horizontal derivatives. A sensitivity analysis then evaluates the accuracy and stability of the RBF-FD method as a function of the RBF parameters for a representative test function (Section 4). These tests help to identify the optimal method parameters for simulating waves. Then, the 3D version of the model is validated by simulating two sets of experiments conducted in wave basins with variable bathymetries (Section 5). Finally, an analysis of the advantages and disadvantages of using the RBF-FD approach to simulate nonlinear wave propagation is summarized, including propositions for future work (Section 6).

Brief review of RBF methods

RBFs were first introduced by [START_REF] Hardy | Multiquadric equation of topography and other irregular surfaces[END_REF] for interpolating surfaces. A function f (x) is sought to approximate a given function f (x), for which the values f i (i = 1, ..., N ) are known for a given set of N data points x i (i = 1, ..., N ) using a set of basis functions ψ i (x) (i = 1, ..., N ):

f (x) ≈ N i=1 λ i ψ i (x), (1)
with x = (x, y). The interpolation coefficients λ i are determined by solving the set of linear equations obtained by enforcing the interpolation conditions: f (x i ) = f i for i = 1 to N . For one-dimensional data, a variety of different basis functions (independent of the data points) lead to non-singular linear systems as long as the N data points are distinct (i.e. Fourier and polynomials series). However, this property is no longer guaranteed when the problem is extended to two dimensions [START_REF] Wright | Radial Basis Function Interpolation : Numerical and Analytical Developments[END_REF]. [START_REF] Hardy | Multiquadric equation of topography and other irregular surfaces[END_REF] proposed using a basis of functions composed of a single radial function ψ(r) centered at each data point :

ψ i (x) = ψ(||x -x i ||) (i = 1, ..., N ).
The initial work of [START_REF] Hardy | Multiquadric equation of topography and other irregular surfaces[END_REF] used the multiquadric (MQ) radial function ψ(r) = √ r 2 + C 2 , with r the radial distance from the center x i to node x and C a strictly positive shape parameter to obtain a continuously differentiable basis function (even when r = 0), where C controls the narrowness of the RBF. A wide variety of radial functions may be used. In a study of scattered data interpolation, [START_REF] Franke | Scattered Data Interpolation: Tests of Some Methods[END_REF] tested 29 interpolation methods for six different test functions. The results showed that MQ functions were among the most accurate, together with inverse multiquadric (IMQ) and thin plate spline (TPS) functions (Table 1).

Overall, the RBF method demonstrated good results for spatial interpolation. Therefore the approach was tested further for estimating derivatives in domains with scattered nodes. [START_REF] Stead | Estimation of gradient from scattered data[END_REF] compared the errors obtained in estimating partial derivatives when using MQ or weighted least square quadratic approximate interpolants. Since RBF interpolants do not have polynomial precision (except when polynomial terms are added), the author recommended using the MQ interpolant for surfaces with significant curvature. Later, [START_REF] Kansa | Multiquadrics -A scattered data approximation scheme with application to computational fluid-dynamics -II Solution to parabolic, hyperbolic and elliptic partial differential equations[END_REF] was the first to use the MQ function to solve partial differential equations (PDEs), namely a Poisson equation, with a collocation method. The coefficients of the RBF approximation of the solution are found by solving the linear system obtained by applying the differential operators to the interpolant for the interior nodes and boundary conditions to the interpolant for the boundary nodes. The resulting matrix is not symmetric and is not proven to be unisolvent. Several improvements to this approach were made, by recovering the symmetry of the matrix using Hermite interpolation to modify the basis functions [START_REF] Fasshauer | Solving partial differential equations by collocation with Radial Basis Functions[END_REF], or by imposing both the PDE and the boundary conditions at boundary nodes to increase the constraints where errors are larger [START_REF] Fedoseyev | Improved multiquadric methods for elliptic partial differential equations via PDE collocation on the boundary[END_REF].

Many different functions can be used in the RBF approach, and as shown in the preliminary study by [START_REF] Franke | Scattered Data Interpolation: Tests of Some Methods[END_REF], the results are sensitive to the choice of function. The most commonly used functions (Table 1) can be broken into two families: piecewise-smooth (PS) and infinitely-smooth (IS) functions.

Name (Acronyme)

Function ψ(r) Condition Regularity Polyharmonic Spline (PHS) r m m odd integer PS Thin Plate Spline (T P S) r m log r m even integer PS For IS-RBFs (with a shape parameter C), the interpolation system will not be singular if the scattered nodes are distinct. PS-RBFs do not depend on a shape parameter, but they present a singularity at the origin. To ensure the unique solvability of the linear system, the interpolant has to be modified by including polynomial terms, requiring additional constraints for the linear system to be well-posed. In this case, the non-singularity of the matrix becomes more restrictive since it requires that the nodes are not only distinct but are also unisolvent in the appended polynomial space (Fornberg and Flyer, 2015a).

Multiquadric (MQ) √ r 2 + C 2 C > 0 IS Inverse Multiquadric (IMQ) 1 √ r 2 +C 2 C > 0 IS Inverse Quadratic (IQ) 1 r 2 +C 2 C > 0 IS Gaussian (GA) e -r 2 /C 2 C > 0 IS
Global vs. local methods. The RBF method was first introduced as a global method, in which estimates at each node depend on all nodes in the domain, leading to a full coefficient matrix. When the matrix becomes too large, it often becomes ill-conditioned. The interpolation coefficients (λ i , Eq.1) become oscillatory with large magnitudes that may lead to a poor interpolation because of numerical cancellations. In this case, the size of the matrix can be reduced by considering smaller domains using domain decomposition algorithms [START_REF] Beatson | Fast Solution of the Radial Basis Function Interpolation Equations: Domain Decomposition Methods[END_REF][START_REF] Wong | Multizone decomposition for simulation of time-dependent problems using the multiquadric scheme[END_REF][START_REF] Zhou | Overlapping domain decomposition method by radial basis functions[END_REF]. In the limit, one can instead use a local approach by defining stencils centered at each node of the domain and including only the N sten -1 nearest neighbors (for a total of N sten nodes in each stencil) in the estimation of the function (Eq.1). [START_REF] Tolstykh | On using radial basis functions in a 'finite difference mode' with applications to elasticity problems[END_REF] were the first to consider this method to estimate derivatives with a RBF-FD approach, followed shortly by [START_REF] Wright | Radial Basis Function Interpolation : Numerical and Analytical Developments[END_REF] and [START_REF] Shu | Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier-Stokes equations[END_REF]. [START_REF] Wright | Scattered node compact finite differencetype formulas generated from radial basis functions[END_REF] improved the accuracy of this method by including a linear combination of derivatives of the function at the surrounding nodes. The local method has the advantage of reducing considerably the computational time in comparison to the global method, as well as being parallelized easily. The construction of approximate formulas for the derivatives using RBF interpolants, also called the RBF-FD method, will be presented in more detail in Section 3.3.

However, even when using local methods, the matrix may become illconditioned for IS-RBFs when C → ∞ (i.e. in the limit of flat basis functions). Several specific algorithms have been developed to obtain accurate results even for large values of C (i.e. Contour Padé (Fornberg and Wright, 2004), RBF-QR [START_REF] Fornberg | A stable algorithm for flat radial basis functions on a sphere[END_REF], RBF-GA [START_REF] Fornberg | Stable calculation of Gaussian-based RBF-FD stencils[END_REF] or more recently the RBF-RA algorithm [START_REF] Wright | Stable computations with flat radial basis functions using vector-valued rational approximations[END_REF]). Using a stable algorithm not only improves the derivative estimation accuracy by allowing the use of larger values of C, but also makes the choice of optimal C less critical. Nevertheless, the algorithms may be difficult to adapt to a specific mathematical model or the modified RBFs may have much more complicated expressions.

PS vs. IS functions. The choice of RBF to obtain the most accurate estimates is not straightforward, and some trade-offs must be considered when using IS or PS functions. Errors when using either family of functions in the RBF method depend on the specific choice of the radial function (e.g. Table 1), the mesh resolution and spacing, and the stencil size N sten . Two main characteristics differentiate these two families of functions: (1) the type of convergence as a function of the node spacing, and (2) the dependency on a shape parameter C.

When using global RBF methods, IS-RBFs have spectral convergence, while PS-RBFs have only algebraic convergence, which often leads to a preference for IS-RBFs (Fornberg and Flyer, 2015a). However, when using local RBF methods (i.e. RBF-FD method), the spectral accuracy of IS-RBFs is lost, minimizing their advantage over PS-RBFs with respect to convergence. Additionally, stagnation errors exist. One type of PS-RBFs, PHS, require the addition of a polynomial (with M terms) to the interpolant to guarantee the unisolvency of the system. The PHS RBFs then have a convergence rate corresponding to the degree of the added polynomial. For complex applications, [START_REF] Barnett | A Robust RBF-FD Formulation based on Polyharmonic Splines and Polynomials[END_REF] and [START_REF] Bayona | On the role of polynomials in RBF-FD approximations: II. Numerical solution of elliptic PDEs[END_REF] recommended an added polynomial of degree such that there are approximately twice as many RBFs as polynomial terms in the interpolant (i.e N sten ≈ 2M ).

With IS-RBFs, the estimation error depends strongly on the value of the shape parameter C: for small values of C, the error is generally high, decreasing with an increase in C. A minimum is often reached for an intermediate value of C = C opt (called the optimal value of C in the following). When C is increased beyond this optimal value, the error increases and large oscillations may be observed as the matrix becomes ill-conditioned. The matrix may already become ill-conditioned for values of C smaller than C opt : in this case, the minimal error is just at the limit of ill-conditioning. In the limit of C → ∞, Fornberg et al. (2004) showed that when the interpolant limit is finite, it tends to a multivariate polynomial. Finding C opt is a difficult task. For global methods, there is no existing mathematical theory to determine the optimal choice of C. Often, this choice is based on the inter-node spacing for convenience. However, based on a series of tests, [START_REF] Carlson | The parameter r 2 in multiquadratic interpolation[END_REF] and then [START_REF] Rippa | An algorithm for selecting a good value for the parameter c in radial basis function interpolation[END_REF] concluded that the optimal value of C depends on the shape of the interpolated function and not on the node positions or spacing. When the RBF method is used to interpolate data, the value of C is chosen by cross-validation methods. For example, [START_REF] Rippa | An algorithm for selecting a good value for the parameter c in radial basis function interpolation[END_REF] developed a method based on the minimization of a cost function calculated as the error between the interpolant and the desired function. [START_REF] Fasshauer | On choosing "optimal" shape parameter for RBF approximation[END_REF] adapted this algorithm for the resolution of PDEs with RBF pseudospectral methods. For RBF-FD, [START_REF] Bayona | RBF-FD formulas and convergence properties[END_REF] derived an expression of the estimated error as a function of C, showing that C opt depends on the shape of the function and its derivatives. They also showed that it is independent of the node spacing at first order, but can vary with node locations in 2D. Given an expression of the estimated error, [START_REF] Bayona | Optimal constant shape parameter for multiquadric based RBF-FD method[END_REF] proposed an algorithm to find the optimal value of C. However, this requires first estimating the derivatives at each point with another less accurate method (since the values of the derivative of the function are necessary to compute the error estimates), before re-estimating the derivatives with the RBF-FD method and newly obtained optimal value of C.

The selection of an optimal or even "good" value for the shape parameter C can be challenging. Thus, RBFs without shape parameters, such as PHS, recently have become more attractive, noting furthermore that they produce relatively well-conditioned matrices [START_REF] Bayona | RBF-FD formulas and convergence properties[END_REF][START_REF] Flyer | On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy[END_REF].

For more details on RBF methods and numerous application examples, see the recently published book of Fornberg and Flyer (2015a).

Model description

Mathematical model

The fluid domain is delimited in the vertical by the free surface at z = η(x, t) and the bottom at z = -h(x), with the vertical axis z pointing upwards and the origin at the still water level. The fluid is accelerated by gravity g. At the free surface, the atmospheric pressure is assumed uniform and constant in time (chosen to be 0 by convention), and surface tension is neglected. Potential flow theory is adopted by further considering the irrotational flow of an inviscid and homogeneous fluid of constant density. The kinematics of the fluid are obtained from the velocity potential Φ(x, z, t) such that v = (u, v, w) = (∇Φ, Φ z ), where ∇f ≡ (f x , f y ) is the horizontal gradient operator, and partial derivatives are denoted with subscripts.

If the free surface is assumed to be single-valued in x (no overturning waves), the evolutions of η(x, t) and Φ(x, t) ≡ Φ(x, z = η(x, t), t) are governed by the nonlinear kinematic and dynamic free surface boundary conditions, expressed as functions of free surface quantities only [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF]:

η t = -∇η. ∇ Φ + w 1 + (∇η) 2 , (2) Φt = -gη - 1 2 (∇ Φ) 2 + 1 2 w2 1 + (∇η) 2 , (3)
where w(x, t) ≡ Φ z | z=η(x,t) is the vertical velocity at the free surface.

To integrate these equations in time, w is determined from η and Φ by solving the Laplace BVP for the velocity potential Φ in the entire domain (Laplace equation) supplemented by the free surface and bottom boundary conditions (BCs):

∇ 2 Φ + Φ zz = 0, -h(x) ≤ z ≤ η(x, t), (4) Φ = Φ(x, t), z = η(x, t), (5) ∇Φ • ∇h + Φ z = 0, z = -h(x). (6)
At the lateral boundaries, periodic, Dirichlet or Neumann BCs are imposed to close the problem. The Laplace BVP is solved using a spectral approach in the vertical dimension following the work of [START_REF] Tian | A numerical model on the interaction between nearshore nonlinear waves and strong currents[END_REF]. The method is described for the case of a single horizontal dimension (i.e. x = x) in [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF] and [START_REF] Raoult | Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments[END_REF]. The extension to two horizontal dimensions is quite straightforward, and only the main steps are presented here.

First, a change of the vertical coordinate from z ∈ [-h(x), η(x, t)] to s ∈ [-1, 1] is made to project the time-varying domain to a constant height domain extending from the bottom at s = -1 to the free surface at s = +1:

s(x ¯, z, t) = 2z + h -(x ¯, t) h + (x ¯, t) , (7)
where

h + (x ¯, t) ≡ h(x ¯) + η(x ¯, t) and h -(x ¯, t) ≡ h(x ¯) -η(x ¯, t).
The BVP is then reformulated in the transformed space (x, s) for Φ(x, z, t) ≡ ϕ(x, s(x, z, t), t).

Second, the vertical variation of the velocity potential is approximated by a linear combination of Chebyshev polynomials of the first kind, T n (s):

ϕ(x, s) ≈ N T n=0 a n (x)T n (s), (8)
where n = 0, 1, 2, ... indicates the order of the polynomial, and N T is the maximum order of the Chebyshev polynomials. These polynomials are easy to compute, form an orthogonal basis over the range [-1, 1], and converge rapidly over a large domain. [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF] and [START_REF] Raoult | Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments[END_REF] have shown that values of N T smaller than 10 (typically in the range [5,[START_REF] Amestoy | A fully synchronous multifrontal solver using distributed dynamic scheduling[END_REF]) are sufficient to obtain high accuracy for a variety of 2DV nonlinear wave propagation test cases. The approximation Eq.( 8) is then introduced in the BVP, and a Chebyshev-Tau method [START_REF] Boyd | Chebyshev and Fourier Spectral Methods: Second Edition, Revised[END_REF] is applied to eliminate the dependence on the vertical coordinate s. The Laplace equation ( 4) is projected on polynomials T p for p ranging from 0 to N T -2, supplemented by the boundary conditions, Eq.( 5) and ( 6), respectively. For nodes located on a lateral boundary, the Laplace equation ( 4) is replaced, either by a Neumann condition for an impermeable boundary ∇Φ.n lat = 0 (where n lat is the unit normal vector at the lateral wall) or by a Dirichlet condition for wave generation, imposing the velocity potential from linear theory (generally supplemented with a relaxation zone). The final set of equations is a linear system for the a n (x), which depends only on x at a given time:

                                               a p,xx + a p,yy + N T n=0 C x pn a n,x + N T n=0 C y pn a n,y + N T n=0 D pn a n = 0 for p = 0, ..N T -2 N T n=0 a n = Φ(x ¯, t) h + h x N T n=0 (-1) n a n,x + h + h y N T n=0 (-1) n a n,y + 2 (1 + h 2 x + h 2 y ) N T n=0 (-1) n-1 n 2 a n = 0 (9a) (9b) (9c)
where

C x pn = (m 0101 B p01n + m 1101 B p11n )/m 0220 , (10)
C y pn = (m 0011 B p01n + m 1011 B p11n )/m 0220 , (11)
D pn = (m 0002 B p02n -m 1002 B p12n + m 2002 B p22n + m 0001 B p01n + m 1001 B p11n )/m 0220 .(12)
The m ijkl terms depend only on h + , h -, and their spatial derivatives (up to order two). The expressions for these terms are shown in Appendix A. The B pikn terms are introduced to express the projection of the terms s i d k Tn ds k on the polynomial T p , and they can be determined analytically as a function of n and p (see Appendix B). These terms are constant and can be computed once at the beginning of each simulation, after N T is chosen.

For each node x, N T + 1 unknown coefficients a n in Eq.( 8) must be determined. With the RBF-FD method, the horizontal spatial derivatives of the a n coefficients are approximated as linear combinations of the values in the vicinity of the node considered, leading to a set of coupled linear equations.

Once the a n (x) coefficients are determined, the vertical velocity at the free surface w(x) is obtained readily from:

w(x) = Φ z (x ¯, z = η) = s z ϕ s (s = +1) ≈ 2 h + (x ¯) N T n=1 a n (x ¯) n 2 (13)
The vertical velocity can then be used to evaluate the right hand side of Eq.( 2) and Eq.( 3), required by the numerical scheme to integrate in time.

Numerical implementation

A classical, explicit fourth-order Runge Kutta (RK4) scheme with a constant time step is used to integrate Eq.( 2) and Eq.(3) in time. At each sub-step of the RK4 scheme, the Laplace BVP is solved. The domain is discretized with N P XY nodes in the horizontal (x, y)-plane, and N T is the maximum order of the Chebyshev polynomials in Eq.( 8). Horizontal derivatives are approximated with the local RBF-FD method (described in section 3.3) using a stencil consisting of a fixed number of neighboring nodes. Therefore, for nodes far from the boundaries, the stencil dimensions are symmetric, while they become asymmetric at and near the boundaries. The effect of this asymmetry will be studied further in section 4. The RBF-FD method is chosen for its ease of implementation, with an algorithm similar to finite difference methods, and for its flexibility with scattered nodes enabling the simulation of complicated domain geometries and the possibility of local refinement. The linear system corresponding to the Laplace BVP is composed of N P XY (N T + 1) equations for the coefficients a n (x i ), with n = 0, ..., N T and i = 1, ..., N P XY . The associated matrix is sparse, and the system is currently solved using the direct solver MUMPS ("MUltifrontal Massively Parallel Solver", v4.10.0) [START_REF] Amestoy | A fully synchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF], using the default settings.

RBF-FD method 3.3.1. Theoretical background

Similar to finite difference methods, applying a linear differential operator L to a given function f at x 1 is expressed as a linear combination of the values of the function f at the nodes in the stencil:

L f (x 1 ) ≈ Nsten i=1 α L i f (x i ). ( 14
)
The stencil is composed of the node itself and its N sten -1 nearest neighbors (x 2 , x 3 , ..., x Nsten )

The weights α L i , for i = 1 to N sten , are determined by requiring the approximation in Eq.( 14) to be satisfied by the set of radial functions centered at each node of the stencil

ψ i (x) ≡ ψ(||x -x i ||), i ∈ [1, N sten ]
. This leads to the resolution of a linear set of N sten equations for N sten unknowns, which can be written as:

     ψ(||x 1 -x 1 ||) • • • ψ(||x Nsten -x 1 ||) ψ(||x 1 -x 2 ||) • • • ψ(||x Nsten -x 2 ||) . . . . . . . . . . . . ψ(||x 1 -x Nsten ||) • • • ψ(||x Nsten -x Nsten ||)         α L 1 . . . α L Nsten    =    L ψ(||x -x 1 ||)(x 1 ) . . . L ψ(||x -x Nsten ||)(x 1 )    (15)
As mentioned in the introduction, several different forms of RBF ψ(r) can be used (Table 1).

To improve the accuracy of the estimation, especially at the boundaries of the domain [START_REF] Fornberg | Observation on the behaviour of radial basis function approximation near boundaries[END_REF], the RBF can be supplemented with a polynomial of degree l: M j=1 b j p j (x), where p j (x) M j=1 is a basis of polynomials up to degree l in R 2 and M = l+2 2 . In this case, M additional constraints minimizing the far-field growth must be enforced to close the system:

Nsten i=1 λ i p j (x i -x 1 ) = 0 j = 1, 2, 3..., M. (16)
Finally, the system to be solved to obtain the α L i coefficients becomes:

Γ P P T 0 α β = L ψ L p ( 17
)
where Γ is the matrix of ψ in the left hand side of Eq.( 15), P is a M × N sten matrix formed by the p j , j = 1, ..., M basis of polynomials up to degree l in R 2 , and β and L p are the M × 1 vectors formed by β i and L p i (x 1 ). The application of the method depends on the non-singularity of the matrix Γ P P T 0 . The matrix Γ is guaranteed non singular for IS-RBFs provided that the nodes are distinct. For PS-RBFs Γ is no longer guaranteed to be nonsingular, and a polynomial of degree l must be added to ensure that the system is uniquely solvable [START_REF] Wright | Radial Basis Function Interpolation : Numerical and Analytical Developments[END_REF]. [START_REF] Barnett | A Robust RBF-FD Formulation based on Polyharmonic Splines and Polynomials[END_REF] showed that for PHS of the form r m , the degree l of the added polynomial must satisfy l ≥ (m -1)/2. The addition of a polynomial also requires the use of an unisolvent set of nodes [START_REF] Wright | Radial Basis Function Interpolation : Numerical and Analytical Developments[END_REF].

Numerical implementation

The weights α L i must be calculated at all nodes in the domain for each differential operator L required for the discretization of the PDEs. In the present model, first and second-order derivatives in the two horizontal dimensions are calculated. At each node x of the horizontal mesh, the following steps are carried out. First, the N sten -1 nearest neighbors of the selected node are identified. The size of the stencil and the degree of the augmented polynomial are defined at the beginning of the simulation and are constant for all nodes. Then, the linear system Eq.( 17) is solved by completing a LU decomposition. The weights are computed once at the beginning of the simulation and are subsequently used to estimate all spatial derivatives in the model, including, but not limited to, the free surface elevation η, the free surface potential Φ, the water depth h, and the a n coefficients.

The implementation of the RBF-FD method is first tested to evaluate its accuracy in approximating spatial derivatives of a representative wave function (Section 4), and then it is applied for the simulation of wave tank experiments with variable bathymetries (Section 5).

Accuracy of the RBF-FD derivative estimates

Method

A series of tests were conducted to evaluate the capabilities and limitations of the RBF-FD method in estimating first and second-order spatial derivatives. These tests evaluate the impact of several parameters such as the RBF type, the value of the shape parameter C in the case of IS-RBFs, the degree of the added polynomial, and the stencil size N sten . Previous work (e.g. [START_REF] Franke | Scattered Data Interpolation: Tests of Some Methods[END_REF]; [START_REF] Stead | Estimation of gradient from scattered data[END_REF]) has tested different types of RBFs (Table 1) to evaluate the interpolation and derivate estimate accuracy. However, these functions were usually chosen arbitrarily, with more or less complex spatial variations. When simulating ocean wave propagation, the free surface generally has oscillatory variations, thus a sinusoidal function is chosen here as a basic representative model:

f (x, y) = A cos 2π L (x cos θ + y sin θ) , (18)
where L = 0.5 m is the characteristic length of variation (or wavelength), A is the wave amplitude such that A/L = 0.05, and θ = 20°is the direction of wave propagation with respect to the x axis. The domain is defined as 0 ≤ x ≤ 1 m and 0 ≤ y ≤ 1 m and is discretized with a regular set of nodes with node spacing ∆x = ∆y = 0.05 m (= L/100). First and second-order derivatives in both horizontal dimensions are estimated with the RBF-FD method and compared to the analytical values (denoted Lf theo hereafter).

The accuracy of the estimation is evaluated by calculating the normalized averaged error for all N nodes in the domain:

RM S Error = N i=1 (Lf (x i ) -Lf theo (x i )) 2 N i=1 (Lf theo (x i )) 2 (19)
Globally, errors are larger closer to the boundaries, and more particularly when the stencil is one-sided in the direction of the derivative (i.e. for x = 0 m and x = +1 m for x-derivatives, and y = 0 m and y = +1 m for y-derivatives). Two subsets of nodes are then defined, based on the asymmetry of the stencil: the interior nodes with a centered stencil (such as node A in Figure 1) and the boundary nodes with asymmetric or one-sided stencils (such as nodes B and C in Figure 1). In the following, the global averaged error Eq.( 19) is shown except when the error behavior differs for the two node sets and is analyzed separately.

From the series of tests carried out to study the sensitivity of the method's accuracy to the chosen parameters, one initial question arises: "Which RBF is optimal?". The values of the other parameters will then depend on this choice. For IS-RBFs, which depend on a shape parameter, the accuracy of the method was evaluated as a function of the: IS-RBF function: MQ, IMQ, IQ and GA (see Table 1), variation of the shape parameter in the range C ∈ [0, 30] m, inclusion of an added polynomial up to degree 2, and stencil size for N sten = 5, 13 and 21, corresponding to the optimal thresholds defined by [START_REF] Bayona | RBF-FD formulas and convergence properties[END_REF] for regular node sets.

For PS-RBFs (not depending on a shape parameter), the accuracy of the method was evaluated as a function of the: PS-RBF function: one T P S function ψ(r) = r 4 log r to have at least the continuity of the second-order derivatives, and two PHS functions ψ(r) = r 5 and ψ(r) = r 7 , inclusion of an added polynomial up to degree 5, (in the following, the shorthand form r 5 + p3 is used to denote the PHS r 5 with an added polynomial of degree 3), and stencil size in the range N sten ∈ [9, 56].

Results

Sensitivity to the shape parameter C

In Figure 2, the global averaged error for the four IS-RBFs is plotted as a function of the shape parameter C for the first-order (f x , f y ) and second-order (f xx , f yy ) derivatives. For IS-RBFs, the derivative estimate accuracy depends strongly on C. The four RBFs display the same general behavior as a function of C: large errors for small values of C that decrease with an increase in C, and the development of oscillations when C exceeds an unknown threshold (depending on the RBF). These oscillations appear when the coefficient matrix (Eq.( 17)) becomes ill-conditioned. The optimal value of C for which the error is minimum (denoted C opt ) depends on the RBF and derivative estimated, and is therefore not known a priori (as previously stated in Section 2). First and second-order derivatives in x reach a minimum for C ≈ 0.42 m, whereas the matrix becomes ill-conditioned for the first and second-order derivatives in y before a minimum is observed. It can be inferred that the optimal value of C for derivatives in y is larger than for derivatives in x for this particular function and value of θ. The choice of the value of C is thus a compromise between optimizing the accuracy of the solution and reaching the limit of an ill-conditioned matrix. The accuracy of the estimation could be improved by allowing the shape parameter to vary with the RBF center (i.e. [START_REF] Kansa | Improved accuracy of multiquadric interpolation using variable shape parameters[END_REF][START_REF] Kansa | Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations[END_REF], which produces larger variations in the matrix coefficients, thus reducing the condition number. [START_REF] Bayona | Optimal variable shape parameter for multiquadric based RBF-FD method[END_REF] developed an algorithm for the RBF-FD method to find the C opt for each node, but it is an inefficient approach in the current model since it requires estimating the derivatives twice. Finally, C opt depends on the function whose derivatives are estimated, suggesting that different coefficients would need to be calculated for each variable. In addition, some variables in the model (η, Φ, ...) are time dependent, requiring the coefficients to be a function of time as well. To increase the accuracy, a simpler approach was tested by normalizing the stencil, or transforming the local support to a unit circle, following the work of [START_REF] Shu | Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier-Stokes equations[END_REF]. This was expected to reduce the difference in accuracy between interior nodes (with centered stencils) and boundary nodes (with asymmetric stencils), allowing the use of larger values of C for nodes with asymmetric stencils. However, tests with normalized stencils did not improve significantly the RBF accuracy (results not shown here).

The errors for the three PS-RBFs do not depend on a shape parameter and are therefore constant as a function of C (Figure 2). For PS-RBFs, the errors decrease when the degree of r increases, and ψ(r) = r 7 produces the smallest errors for the four derivatives considered. Overall, PS-RBFs cause smaller errors than IS-RBFs for small values of C. The RBF ψ(r) = r 7 produces errors comparable to the minimum errors obtained with the IS-RBFs for f y and f yy , and slightly larger errors than what can be reached by the IS-RBFs near C opt . Although the RBF ψ(r) = r 7 may not be the optimal choice for the presented derivatives, it offers the advantage of not relying on the choice of a shape parameter, whose selection is not a straightforward process, as previously mentioned. Finally, since the errors with PHS r 5 and r 7 are smaller than those with r 4 log(r), only the PHS RBFs will be examined in the following tests.

Sensitivity to the stencil size and added polynomial degree

The dependence of the estimation accuracy on the stencil size (N sten ) and the degree of the added polynomial (l) demonstrates that the sensitivity is not the same for IS and PS-RBFs.

IS-RBFs. First, focusing on IS-RBFs, Figure 3 shows the global error for the estimation of f x as a function of C for three values of the stencil size, N sten = 5, 13, and 21. The accuracy of the derivative estimates improves significantly by increasing the stencil size from 5 to 13 nodes, and even more by increasing to 21 nodes. However, C must be chosen carefully since the range of values of C producing a well-conditioned matrix is reduced. The matrix becomes ill-conditioned for smaller values of C when N sten increases (e.g. C ≈ 0.8 m for N sten = 13, whereas C ≈ 0.5 m for N sten = 21). In addition, increasing N sten increases the computational time(tests of regular wave propagation, not shown here, exhibited a computational time proportional to N 1.36 sten , using a fixed value of N P XY ), so a compromise must be made between the desired accuracy, the difficulties in finding an optimal value of C, and the computational time.

The sensitivity of the error to the degree of the added polynomial is then studied by increasing l from 0 to 2 and comparing the results to those without an added polynomial. The boundary and interior nodes present different behavior as a function of l (Figure 4). Despite the smaller number of boundary nodes (in comparison with interior nodes), the global error (calculated for all nodes, Figure 4c) is dominated by the boundary errors, therefore presenting the same dependence on the shape parameter C. For interior nodes (Figure 4a), adding higher degree polynomials decreases the error for C < 0.18 m, but this improvement is lost for higher values of C. For boundary nodes (Figure 4b), the errors are reduced significantly with the addition of a polynomial and an increase in the degree of the polynomial for C < 0.3 m. Since RBFs are not exact approximations of polynomials, it is essential to add at least a constant to the RBF interpolant to be able to estimate accurately the derivative of a constant function.

For IS-RBFs, increasing the stencil size improves the accuracy but causes the matrix to become ill-conditioned for smaller values of C. The addition of high-order polynomials is only beneficial for small values of C, especially for the boundary nodes, in agreement with [START_REF] Fornberg | Observation on the behaviour of radial basis function approximation near boundaries[END_REF]. PS-RBFs. The addition of a high-order polynomial is essential to guarantee the inversibility of the collocation matrix for PHS RBFs. A minimum polynomial degree l is required, which depends on the degree of the PHS function. In addition, the maximum degree of the added polynomial is limited by the size of the stencil. To ensure that the problem is well-posed, N sten must be larger than the number of independent monomials constituting the basis of polynomials of the same degree as the added polynomial (i.e. [START_REF] Barnett | A Robust RBF-FD Formulation based on Polyharmonic Splines and Polynomials[END_REF]. Given these constraints, a series of tests were conducted to study the sensitivity of the error estimation of the PHS r 5 and r 7 to the stencil size (N sten ∈ [9, 56]) and to the degree of the added polynomial (l ∈ [2, 5]).

With regular node sets, the condition on the minimum stencil size for a given degree of added polynomial is not sufficient to ensure the nonsingularity of the matrix. The regularity of the node set does not allow the matrix to be unisolvent for the polynomial basis [START_REF] Wright | Radial Basis Function Interpolation : Numerical and Analytical Developments[END_REF]. The stencil size thus has to be increased to recover the inversibility of the matrix [START_REF] Barnett | A Robust RBF-FD Formulation based on Polyharmonic Splines and Polynomials[END_REF]. With an irregular node set, this may not occur. The results obtained with both PHS (r 5 and r 7 ) are compared in Figure 5 for the first derivative in x. Similar behavior is observed for the two PHS, but PHS r 7 produces smaller errors for a given degree of the added polynomial. The interior node and boundary node errors vary differently as a function of N sten but in both cases, the errors decrease with an increase in the degree of the added polynomial. For an added polynomial with an even degree, the error for the interior nodes is weakly dependent on the stencil size N sten , whereas for an added polynomial with an odd degree, a minimum is obtained for any stencil size with a symmetric distribution of nodes. This effect is caused by the regular distribution of the nodes and may not be observed with irregularly spaced nodes (for example, for boundary nodes with asymmetric stencils). The same trends are observed for the estimation of second-order derivatives, but for the inverse of odd and even added polynomials (not shown here). Contrary to the PHS r 7 , the PHS r 5 has the advantage of being used with only a second degree added polynomial, thus requiring a smaller N sten and allowing a reduction in the computational time.

Convergence study as a function of the grid spacing

For IS-RBFs (e.g. IMQ, Figure 6) the optimal C is generally insensitive to the node spacing (here C opt 0.4 m). As the node spacing decreases, the errors decrease until a certain limit below which the matrix becomes illconditioned. Thus C is generally decreased with the node spacing to keep the condition number of the collocation matrix roughly constant, at about 10 15 according to [START_REF] Flyer | On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy[END_REF]. This causes error saturation and the loss of convergence, which can be restored by adding polynomials.

For PHS RBFs, there are no saturation errors, and the convergence rate depends on the degree of the added polynomial l (∝ ∆x l-k+1 ), where k is the order of the estimated derivative [START_REF] Flyer | On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy[END_REF]. When the degree of the PHS RBF increases (from r 5 to r 7 ), the errors decrease slightly only for low order added polynomials (i.e. l = 3 and 4). The convergence rate is independent of the stencil size N sten (Figure 7).

These results, concerning error estimation as a function of node refinement for IS-RBFs and PS-RBFs (PHS), are in agreement with the studies of [START_REF] Bayona | RBF-FD formulas and convergence properties[END_REF] and [START_REF] Flyer | On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy[END_REF].

Summary of the derivative estimate tests

Finally, based on a literature review and a series of tests including, but not limited to those shown here, the use of the PHS r 7 with an added polynomial of degree 3 and N sten = 21, appears to be a good compromise for applica-tions similar to those presented here. For general applications, determining the optimal shape parameter may be challenging (where the derivatives of all variables are estimated with the same coefficients). It is more practical and efficient to use a RBF that does not depend on a shape parameter since C opt depends strongly on the estimated derivative and the function itself. The PHS r 7 is preferred to the PHS r 5 because, for the same degree of added polynomial, the PHS r 7 tends to produce smaller errors. These choices are recommended to minimize the expected errors and stencil size (i.e. computational time) for the estimation of first and second-order derivatives, and will be used in the subsequent test case simulations. This recommendation is derived from a sensitivity study testing a single function. Although there is no formal proof that these results can be extended to a wide range of functions, the chosen sinusoidal function is assumed representative of the type of functions encountered in real applications, suggesting the generalization of the conclusions concerning the efficiency and accuracy of this method for wave propagation simulation models.

Validation test cases

The 3D version of the model is validated by comparing simulation results to measurements from two laboratory experiments studying the convergence of regular and irregular waves propagating over two different bathymetric profiles: a semi-circular step based on the experiments of [START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF], and a submerged shoal based on the experiments of [START_REF] Vincent | Refraction-diffraction of irregular wave over a mound[END_REF].

Nonlinear wave propagation over a semi-circular step

Whalin (1971) performed a series of experiments of regular waves propagating over a semi-circular bottom topography that acts as a focusing lens. These experiments were conducted to test the limit of linear and nondiffractive theory in a convergence zone, considering non-breaking waves with periods of 1, 2, and 3 s, for three wave heights. The bottom topography was designed to produce strong wave convergence along the centerline of the basin and to minimize sidewall effects and dissipation by bottom friction. The wave tank was 6.096 m wide and 25.603 m long. In the experiments, regular waves were generated with a piston wave maker and propagated from an initial water depth of h 0 = 0.4572 m to a shallower region of depth h 1 = 0.1524 m. The bathymetric profile (Figure 8) is defined analytically by:

h(x, y) =      h 0 , -20.0 ≤ x ≤ 10.67 -G(y) h 0 + 1 25 (10.67 -G(y) -x), 10.67 -G(y) < x < 18.29 -G(y) h 1 , 18.29 -G(y) ≤ x ≤ 35 (20)
with G(y) = y(6.096y), (h, x and y in m). Four sets of wave conditions were simulated (see Table 2). For wave condition A, the computational domain extends from -3.91 m to 32.5 m in the x-direction and from 0 to 6.096 m in the y-direction. The computational domain is longer than the physical domain in the x-direction to include a one-wavelength long relaxation zone at the left boundary for wave generation and for the absorption of waves reflected from the underwater topography. A three-wavelength long relaxation zone is added at the right boundary for wave absorption. Impermeable conditions are applied at the lateral boundaries.

The domain is discretized with N P XY = 137, 712 scattered nodes with regular node spacing (∆x ≈ ∆y ≈ 0.04 m, or approximately L/98). Waves are propagated during 36 s (18 periods), with a constant time step ∆t = 0.0178 s (≈ T /112) and N T = 7. [START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF], where the 0 subscript denotes deep water conditions.

Wave condition T (s) A (m) L 0 (m) k 0 h 0 k 0 A A 2 0.
The free surface profile at the end of the simulation (Figure 9) shows a quasi-2D behavior with almost no variations in the y-direction in the deeper part of the domain (x < 7.5 m). The 3D wave patterns develop in the shallower zone where nonlinear effects are important. The convergence of wave energy is caused by a combination of shoaling, diffraction, and refraction over the convergent bathymetric profile. The simulated crest and trough elevation envelopes along the centerline of the tank (y = 3.048 m) show that before the foot of the slope, the crest and trough are nearly symmetric with respect to the still water level (Figure 10). In the shallower zone (x > 15 m), the waves are narrower and their amplitudes increase, with deeper troughs and higher crests, breaking the horizontal symmetry observed in the deeper part of the domain. The vertical asymmetry of the waves also increases around x = 10 m, displaying a steeper wave front. At the maximum of the crest envelope (x ≈ 20 m), the wave has two small lobes on each side, a consequence of the increase of the second harmonic amplitude caused by nonlinear effects on the slope. [START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF]. The wave envelope indicates the maximum and minimum free surface elevation during the simulation (gray lines). The light gray shaded areas (x < 0 m and x > 25 m) indicate the relaxation zones for wave generation and absorption in the numerical model.

To examine more closely nonlinear effects and energy transfers between harmonics, a Fourier analysis of the simulated wave signal was completed along the centerline of the wave tank (Figure 11). The model accurately reproduces the spatial evolution of the amplitudes of the first three harmonics (frequencies f , 2f and 3f ) in comparison to the measurements. The amplitude of the second harmonic is slightly underestimated in the deeper part of the domain (for x ≤ 12 m), which may be related to the linear method used to generate waves in the model. Conversely, in the shallower zone, the second harmonic amplitude is slightly overestimated. As mentioned previously, in the convergence region (around x = 20 m), the second and third harmonic amplitudes increase due to energy transfers from the first harmonic, and the second harmonic amplitude becomes nearly half of the first harmonic amplitude. Despite the energy transfers to higher order harmonics, the amplitude of the first harmonic does not decrease as one would expect in the case of an alongshore uniform bathymetric profile. According to [START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF], along the centerline of the tank, the rate of decrease in amplitude due to nonlinear transfers to higher harmonics is compensated by the rate of increase in amplitude due to refraction and shoaling.

To test the sensitivity and flexibility of the model to the computational domain node distribution, the simulation was also run using two irregular node sets, one with a homogeneous node spacing (≈ 0.06 m) and the second with node spacing varying with the water depth (between ≈ 0.06 m in the shallower part and ≈ 0.1 m in the deeper part). The first irregular node distribution was created with the 3D finite element generator Gmsh (gmsh.info) distributed under GNU GPL. A zoom of a small part of the domain is shown (Figure 12) for the irregular (left) and regular (right) node sets. The irregular node set distribution is not quasi-uniform, however there were no instability problems. For the second irregular node distribution, another node generator (Fornberg and Flyer, 2015b) enabling easy refinement with the bathymetry and maintaining high regularity at local level was used. A repel algorithm, based on [START_REF] Bayona | On the role of polynomials in RBF-FD approximations: II. Numerical solution of elliptic PDEs[END_REF], was also applied after the addition of nodes close to the boundary, in the deeper part, to ensure simulation stability. The resulting harmonic amplitudes are almost superimposed with those obtained with the regular node set for both irregular node sets (Figure 11).

The sensitivity to the spatial resolution was studied by running seven simulations with regular node spacing (∆x = 0.04, 0.06, 0.075, 0.09, 0.16, 0.24, 0.32 m). The resulting harmonic amplitudes are nearly indistinguishable for ∆x ≤ 0.16 m (Figure 13) with normalized root mean square errors less that 1.5%, 4% and 7% for the first three harmonics, respectively. For ∆x > 0.16 m, larger differences appear for all three harmonics, in particular in the shallow water region over the step (x > 15 m). For simulations with ∆x ≥ 0.09 m, instabilities sometimes appeared at the boundaries. Therefore, in these simulations, the internodal distance was decreased near the boundary to ensure stability. To evaluate the computational efficiency as a function of the number of nodes discretizing the computational domain (N P XY ), the computation time of one iteration is shown for each simula- tion (Figure 14). The CPU time is proportional to N P XY 1.76 . This is lower than the cost of the traditional exact sparse factorization (N P XY 2 ), which may be due to the use of a multifrontal factorization method in the MUMPS solver. This demonstrates the advantages of minimizing N P XY (Figure 13) while maintaining satisfactory results. [START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF].

The simulated spatial evolution of the first three harmonics along the centerline of the tank are compared to the experimental data for the other three wave conditions (Table 2, Figure 15). For wave condition B (Figure 15a), corresponding to the same wave period as case A but with a larger wave amplitude, nonlinear effects are more important, and the second harmonic amplitude is almost two-thirds of the maximum of the first harmonic amplitude. The amplitude of the first harmonic also decreases slightly around x = 20 m, suggesting that the nonlinear energy transfers to higher frequencies occur at a faster rate than the energy convergence from refraction [START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF]. The results of the simulation of wave condition D (Figure 15c), which is less nonlinear but has more important dispersive effects, are in close agreement with the experiments. The last test, wave condition C (Figure 15b), corresponds to nearly shallow water conditions. The simulation overestimates the first harmonic amplitude, but underestimates the second and third harmonic amplitudes. This behavior has been observed in previous studies using a variety of different numerical models [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[END_REF]; [START_REF] Beji | A formal derivation and numerical modelling of the improved Boussinesq equations for varying depth[END_REF]; [START_REF] Engsig-Karup | An efficient flexibleorder model for 3D nonlinear water waves[END_REF]; [START_REF] Wu | Efficient nonhydrostatic modeling of surface waves from deep to shallow water[END_REF]; [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF]; [START_REF] Filippini | Free surface flow simulation in estuarine and coastal environments: numerical development and application on unstructured meshes[END_REF]). [START_REF] Kazolea | An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations[END_REF] suggested that the differences may be caused by the propagation of free reflected waves in the tank and the increased complexity of the case due to the shorter evolution distance and the combination of refraction-diffraction and nonlinearities. For case C, the model's sensitivity to the initial wave amplitude or the position and length of the absorption relaxation zone was not able to explain the differences between the experimental data and the numerical results.

Wave refraction and diffraction over an elliptical shoal

The last test case simulates the propagation of regular and irregular waves over a submerged elliptical mound, reproducing the experiments of [START_REF] Vincent | Refraction-diffraction of irregular wave over a mound[END_REF]. The aim of the experiments was to investigate the limits of approximating irregular wave conditions with monochromatic waves, and they produced a large experimental data set of both monochromatic waves and irregular waves with narrow or broad frequency and directional spread- 2) of the experiments of [START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF].

ing. Two cases are considered here : first the regular wave case M 1, with T = 1.3 s, L = 2.3 m, and A = 0.0275 m, and then the irregular wave case U 3, generated with a JONSWAP spectrum with H s = 0.0254 m, T p = 1.3 s, and a peak enhancement factor γ = 2.

The experiments were conducted in a directional wave basin that was 35 m wide and 29 m long. The measurement area was restricted to a 6.10 m wide by 15.24 m long zone around the elliptical shoal, which had a major axis of 3.96 m and a minor axis of 3.05 m, with the center at (x 0 ,y 0 )=(6.10 m, 13.72 m). The shoal boundary (Figure 16) is defined by: S(x, y) =

xx 0 3.05

2 + y -y 0 3.96 2 = 1. ( 21
)
The water depth around the shoal is constant h(x, y) = 0.4572 m (i.e. for S(x, y) > 1), and the water depth over the shoal is:

h(x, y) = 0.9144 -0.7620 1 - x -x 0 3.81 2 - y -y 0 4.95 2 . ( 22
)
The minimum water depth above the center of the shoal is therefore h min = 0.1524 m. In the experiments, waves were generated with a directional wave generator, located at x = 0 m. The free surface elevation was measured using an array of nine parallel resistive probes placed along nine different transects (five parallel and four perpendicular to the wave maker) during nine different experimental runs. In the following, two transects will be studied (Figure 16): the transversal transect 4 (x = 12.2 m) and the longitudinal transect along the centerline of the tank, consisting of transects 7 and 9 (y = 13.72 m).

To limit the computational time, the simulated domain is smaller than the experimental wave basin. The numerical domain extends from -2.3 m ≤ x ≤ 20.5 m and 3.7 m ≤ y ≤ 23.7 m. Two relaxation zones are added (hatched zones in Figure 16): a one-wavelength long wave generation zone at the left side of the domain, and a two-wavelength long absorption zone at the right side of the domain. Impermeable conditions are applied at the lateral boundaries. The domain is discretized with regularly spaced nodes with ∆x = ∆y = 0.075 m, for a total of N P XY = 81, 435 nodes.

For the regular wave case M 1, waves were generated with an amplitude of A = 0.02325 m, using linear wave theory. This value is slightly smaller than the one prescribed to the wave maker in the experiments, but an adjustment of the incident wave amplitude was necessary to obtain a comparable average wave height (H = 0.0445 m) at a reference probe located in an unperturbed zone of the domain upstream of the shoal (x = 3.05 m, y = 21.34 m). Waves are propagated during approximately 100 s (≈ 78T ), with a constant time step ∆t = 0.036 s (≈ T /36), using N T = 5.

The contour plot of the free surface elevation at the end of the simulation, when the periodic steady state is reached, shows the wave pattern that developed around the shoal (Figure 17). The wave height increases behind the shoal (x > 6 m), with complex 2D patterns with strong variations extending in both horizontal directions. The convergence zone along the centerline of the tank is surrounded by rectilinear zones of almost zero amplitude, look- ing like a wake. In addition, the crests and troughs in the y-direction are modulated with a characteristic length scale of approximately 3 m due to reflections from the lateral walls. The use of a computational domain smaller than the experimental one increases the importance of lateral reflections and possibly overestimates this effect.

To compare the simulation results to the experimental data, a zero upcrossing analysis of the free surface elevation time series is completed to compute the average wave height (H m ) along each transect. To conduct the analysis in the same way as for the experiments, a 28-period window of the free surface elevation time series is considered (once steady state is reached, from t = 60 to 96.4 s). Wave height profiles along the perpendicular transects show good agreement with the experimental data (transect 4, Figure 18a). The wave height profile presents a maximum at the center, corresponding to for case M 1 of [START_REF] Vincent | Refraction-diffraction of irregular wave over a mound[END_REF]. The dotted line indicates the limit of the elliptical shoal on the basin's floor.

the center of the shoal (y = 13.72 m), which is more than twice the incident wave height (ratio ≈ 2.03), but is slightly underestimated in the simulations. Moving symmetrically away from the center, two minima are reached, with wave heights less than half the incident wave height (ratios ≈ 0.21 and 0.43, respectively). Farther from the shoal, the wave height is nearly equal to the incident wave height. The simulated wave height profile in the wave propagation direction also agrees well with the experimental measurements (transect 7-9, Figure 18b), with differences slightly larger than those observed along transect 4. In particular, the increase in the wave height between x = 4-6 m, and the small local peak around x = 7.5 m are not reproduced by the numerical model. After the shoal (x = 9 m), the simulated wave height profile shows small oscillations that may be caused by reflections from the relaxation zone that is not perfectly absorbing. The experimental measurements vary between different runs. At a data measurement point (x = 12.2 m, y = 13.72 m) in two transects, the observed wave height is 0.0975 m and 0.104 m during 2 different runs (along transect 4 and transect 7-9, respectively), which is a difference of approximately 6.25%. Although the variability in the measurements at this location cannot be extended directly to the other measurement points, it can be used to estimate the order of magnitude of the experimental errors and variability.

A harmonic analysis was also performed on the simulated free surface time series along transect 7-9 to show the evolution of the first three harmonic amplitudes (Figure 19). Before the shoal (x < 4 m), the waves are only weakly nonlinear, and the second and third harmonic amplitudes increase over the shoal. Over the shoal (x = 6.10 m), the second harmonic amplitude is more than half first harmonic amplitude. This effect is likely caused by wave convergence and the narrowing of the crest over the shoal (Figure 17). After the shoal, the amplitude of the first harmonic more than doubles due to wave energy convergence induced by refraction. Nonlinearities are not significant after the shoal, although a modulation of the amplitude of the second harmonic is clearly observed. Comparison of the simulated and experimental harmonic amplitude evolution was not possible because of the lack of availability of the experimental time series.

For the irregular wave condition case U 3, the computational domain and numerical parameters remained unchanged, except that the maximum order of the Chebyshev polynomial was increased to N T = 7. N T was increased to account for the presence of higher frequencies in the wave spectrum and the need to resolve accurately the dispersion relation for these frequencies. Waves were propagated during approximately 200 s (≈ 153 T p ). The high frequencies propagate at a lower celerity than the peak frequency, hence the transient period is longer than for case M 1.

The contour plot of the free surface elevation at the end of the simulation (t ≈ 153 T p ) for case U 3 shows more complex and irregular 2D patterns than the regular wave case (Figure 20), and the convergence zone is less well-defined. The effects of reflections from the lateral walls are still visible.

To compare to the experimental data (only available for transect 4), the significant wave height is computed from the free surface elevation time series by calculating H s = 4σ η , where σ η is the standard deviation of the free surface elevation. Time series of around one hundred wave periods are used for this analysis. The simulated H s agrees well with the experimental observations (Figure 21). The pattern of the significant wave height profile is similar to the one obtained for the regular wave case, with a maximum at the center due to wave convergence induced by the bathymetry (ratio ≈ 1.97).

Discussion and conclusions

A meshless approach, based on the RBF-FD method, was chosen for the extension of the 2DV version of a highly nonlinear and dispersive potential wave model to 3D domains. This method has the advantage of being both similar to finite difference methods and simple to implement, not requiring major adaptations to the structure of the 2DV code.

A series of sensitivity tests of the RBF-FD parameters were conducted to examine the robustness of this approach for estimating derivatives of a sinusoidal function. This study demonstrates that accurate results can be obtained with an IS-RBF without significant differences between the tested RBFs: MQ, GA, IMQ and IQ (Table 1). However, IS-RBFs depend on a shape parameter controlling the accuracy of the approximation, and the optimal value of this shape parameter depends on the estimated derivative and the form of the function. To avoid these limitations, it is recommended to use PS-RBFs, which do not depend on a shape parameter. The choice of Figure 20: Contour plot of the free surface elevation at the end of the simulation (t ≈ 153T p ) for case U 3 of [START_REF] Vincent | Refraction-diffraction of irregular wave over a mound[END_REF]. The dotted line indicates the limit of the elliptical shoal on the basin's floor.

the degree of the PS-RBF, the size of the stencil and the degree of the added polynomial is a compromise between the accuracy and computational time.

For a targeted stencil size between 20 and 30 nodes, a PHS function of the form r 7 + p3 is recommended.

The application of the 3D version of the model to simulate two different wave basin experiments showed that complex free surface wave patterns induced by variable bathymetric profiles and the associated nonlinear effects are well reproduced by the model. The nonlinear and dispersive capabilities of the 3D model validate the use of the RBF-FD method for wave propagation in two horizontal directions. Using the RBF-FD approach allows significant flexibility enabling the use of non-rectangular grids and local node refinement, which is of particular interest for applications to real coastal domains. Simulations of wave condition A of [START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF] experiments for two irregular node distributions show that a refinement following the water depth allows increasing node spacing in the deeper part of the domain, thus reducing the number of nodes (from 60,716 to 41,983) and the computational time by approximately 30% while maintaining the difference between harmonic amplitudes below about 2.5%.

With the long term objective of applying the model to real and complex nearshore domains, including wave propagation near coastal and harbor structures, work remains to be done to improve the robustness of the method and the computational efficiency of the numerical model. In some cases, instabilities may occur at or near the boundaries because of derivative estimate errors induced by asymmetric stencils at these locations. In the test cases presented here, refining the mesh close to the boundary was sufficient to avoid the appearance of instabilities. Another option is to add ghost nodes outside of the domain to reduce the one-sidedness of boundary node stencils [START_REF] Barnett | A Robust RBF-FD Formulation based on Polyharmonic Splines and Polynomials[END_REF], but the implementation of such a method is not trivial. More recently, [START_REF] Bayona | On the role of polynomials in RBF-FD approximations: II. Numerical solution of elliptic PDEs[END_REF] showed that the development of instabilities at the boundaries could also be reduced by increasing the stencil beyond the threshold: N sten ≥ 2M . One last method used to stabilize the resolution of PDEs without a physical dissipative term (as is the case of the Zakharov equations) is to add a hyper-viscosity operator to the right hand-side of the evolution equation to introduce artificially a small amount of dissipation that will dampen spurious high frequency oscillations [START_REF] Fornberg | Stabilisation of RBF-generated finite difference method for convective PDEs[END_REF].

The extension of the 2DV code to 3D was accompanied by a significant increase in the computational time. The numerical efficiency of the 3D version of the model needs to be improved in order to perform simulations with large spatial domains. The resolution of the Laplace BVP linear system is the most computationally expensive part of the model, so a parallel version of the linear solver is currently used. Further work could include parallelizing the code with a domain decomposition approach. Other possibilities to reduce the computational cost are currently being explored, including using time integration schemes requiring fewer resolutions of the Laplace BVP (i.e. multi-step predictor-corrector schemes) and/or iterative solvers with suitable preconditioners.

Finally, ongoing work also includes the representation of physical processes, including depth-induced wave breaking and run-up processes to enable simulating nearshore wave environments. Further validation of the 3D version of the model is required, including cases with complex coastlines and variable bathymetries, islands, and coastal structures.

Acknowledgments

Cécile Raoult's PhD thesis was funded partially by the French ANRT (Association Nationale de la Recherche et de la Technologie) with CIFRE grant number 2013-1024.

Figure 1 :

 1 Figure 1: Example of three stencils of N sten = 21 nodes, for a node A in the interior of the domain (symmetric stencil), and nodes B and C on or close to the boundaries (asymmetric stencil).

Figure 2 :

 2 Figure 2: Global errors of the four derivative estimates f x , f y , f xx , and f yy of the test function Eq.(18) for the four IS-RBFs and the three PS-RBFs as a function of the shape parameter C (with N sten = 21 and an added polynomial of degree 0). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Figure 3 :

 3 Figure3: Global error for the estimate of f x for the IMQ RBF as a function of the shape parameter C (with an added polynomial of degree 0) for stencil sizes N sten = 5, 13 and 21.

Figure 4 :

 4 Figure 4: Errors for the estimate of f x , for the IMQ RBF as a function of the shape parameter C (with N sten = 21) for added polynomials of degree 0 to 2: (a) interior nodes, (b): boundary nodes, and (c) all nodes.

Figure 5 :Figure 6 :

 56 Figure 5: Errors for the estimate of f x for PHS r 5 and r 7 as a function of the stencil size N sten and the degree of the added polynomial (shown in the legend) for: (a) interior nodes, (b) boundary nodes, and (c) all nodes.

  p3 Nsten = 20 r 5 + p3 Nsten = 35 r 5 + p3 Nsten = 55 r 7 + p3 Nsten = 20 r 7 + p3 Nsten = 35 r 7 + p3 Nsten = 55

Figure 7 :

 7 Figure7: Errors for the estimate of f x for P HS r 5 and r 7 as a function of the grid spacing (same as in Figure6) for (left) different degrees of the added polynomial l, and (right) several stencil sizes N sten .

Figure 8 :

 8 Figure 8: Bathymetry of the experiments of Whalin (1971).

Figure 9 :

 9 Figure 9: 3D view of the free surface elevation at the end of the simulation (t = 18T ), obtained with PHS r 7 + p3 and N sten = 21 for wave condition A (T = 2 s, A = 0.0075 m and L 0 = 3.91 m) of the experiments of Whalin (1971).

Figure 10 :

 10 Figure10: Free surface elevation profile along the centerline of the tank at the end of the simulation t = 18T (black line) for wave condition A (T = 2 s, A = 0.0075 m and L 0 = 3.91 m) of the experiments of[START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF]. The wave envelope indicates the maximum and minimum free surface elevation during the simulation (gray lines). The light gray shaded areas (x < 0 m and x > 25 m) indicate the relaxation zones for wave generation and absorption in the numerical model.

Figure 11 :

 11 Figure11: Comparison of measured (circles) and simulated spatial evolution of the amplitude of first three harmonics with a regular (solid lines), and two irregular (dashed and dotted lines) node distributions (at frequencies f , 2f , and 3f ) of the free surface elevation for wave condition A (T = 2 s, A = 0.0075 m and L 0 = 3.91 m) of the experiments of[START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF].

Figure 12 :

 12 Figure 12: Zooms of the bottom left corner of the domain for (left) irregular (≈ 0.06 m) and (right) regular (≈ 0.04 m) node distributions used to discretize the computational domain of the experiments of Whalin (1971).

Figure 13 :

 13 Figure13: Comparison of measured (triangles) and simulated (lines) spatial evolution of the amplitude of first three harmonics of the free surface elevation (at frequencies f , 2f , and 3f ) for regular node distributions with different spatial resolution (∆x = 0.04, 0.06, 0.075, 0.09, 0.16, 0.24, 0.32 m) for wave condition A (T = 2 s, A = 0.0075 m and L 0 = 3.91 m) of the experiments of[START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF].

Figure 14 :

 14 Figure 14: CPU time per timestep as a function of the number of nodes discretizing the domain for wave condition A (T = 2 s, A = 0.0075 m and L 0 = 3.91 m) of the experiments of Whalin (1971).

Figure 15 :

 15 Figure 15: Comparison of the observed (circles) and simulated (solid lines) spatial evolution of the amplitude of first three harmonics (at frequencies f , 2f , and 3f ) of the free surface elevation for wave conditions B, C and D (see Table2) of the experiments of[START_REF] Whalin | The limit of applicability of linear wave refraction theory in a convergence zone[END_REF].

Figure 16 :

 16 Figure 16: The numerical domain and bathymetry of the experiments of Vincent and Briggs (1989), with horizontal and vertical white lines indicating the transects along which the simulation results and experimental measurements are compared. Hatched zones indicate the wave generation (left) and absorption (right) zones used in the numerical model.

Figure 17 :

 17 Figure17: Contour plot of the free surface elevation at the end of the simulation (t ≈ 78T ) for case M 1 of[START_REF] Vincent | Refraction-diffraction of irregular wave over a mound[END_REF]. The dotted line indicates the limit of the elliptical shoal on the basin's floor.

Figure 18 :

 18 Figure 18: Average wave height along (a) transect 4 and (b) transect 7-9, for case M 1 of[START_REF] Vincent | Refraction-diffraction of irregular wave over a mound[END_REF]. The horizontal dashed line indicates the incident wave height.

Figure 19 :

 19 Figure19: Simulated spatial evolution of the first three harmonic amplitudes for case M 1 of the[START_REF] Vincent | Refraction-diffraction of irregular wave over a mound[END_REF] experiments along transect 7-9 (y = 13.72 m).

Figure 21 :

 21 Figure 21: Significant wave height along transect 4, for the irregular wave case U 3 of the Vincent and Briggs (1989) experiments. The horizontal dashed line indicates the incident significant wave height.

Table 1 :

 1 Summary of commonly used RBFs describing the function, the constraints on the free parameters, and the regularity of the function (PS stands for Piecewise-Smooth and IS for Infinitely-Smooth).

Table 2 :

 2 Wave characteristics for the four simulations of the experiments of

			0075	3.91	0.7347 0.01205
	B	2	0.0108	3.91	0.7347 0.01736
	C	3	0.0068	6.14	0.4663 0.006936
	D	1	0.0195	1.50	1.9157 0.08171

Appendices

Appendix A m ijkl terms in Eq.( 10), ( 11) and ( 12)

The m ijkl terms appearing in Eq.( 10), ( 11) and ( 12) are terms that only depend on h + = h + η and h -= hη and their spatial derivatives:

Appendix B B pikn terms in Eq.( 10), ( 11) and ( 12)

In Eq.( 10), ( 11) and ( 12) the notation B pikn is introduced:

where f (s) p is the inner product of any function f defined on the interval [-1,1] with a Chebyshev polynomial T p of order p such that :

The following terms have to be estimated: T n p , T n,s p , T n,ss p , s T n,s p , s T n,ss p and s 2 T n,ss p . They can be determined analytically as a function of n and p using the recurrence relation of the Chebyshev polynomials or from linear combinations of previously defined B pikn . Only the final formulas are given below: