# Propagation of 3D nonlinear waves over an elliptical mound with a High-Order Spectral method

Maïté Gouin, Guillaume Ducrozet, P Ferrant

► To cite this version:

Maïté Gouin, Guillaume Ducrozet, P Ferrant. Propagation of 3D nonlinear waves over an elliptical mound with a High-Order Spectral method. European Journal of Mechanics - B/Fluids, 2017, 63, pp.9-24. 10.1016/j.euromechflu.2017.01.002. hal-01437324

HAL Id: hal-01437324

https://hal.science/hal-01437324v1

Submitted on 27 May 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Development and validation of a non-linear spectral model for water waves over variable depth

M. Gouin<sup>a,b</sup>, G. Ducrozet<sup>a</sup>, P. Ferrant<sup>a</sup>

<sup>a</sup>Ecole Centrale Nantes, LHEEA Lab. (ECN, CNRS), Nantes, France

<sup>b</sup>Institut de Recherche Technologique Jules Verne, Bouguenais, France

## Abstract

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.

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

Email address: maite.gouin@ec-nantes.fr (M. Gouin)
# 1. 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 Westphalen et al. (2014), 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 Grilli et al. (2001)) 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 Acquisition technique (see Korsmeyer et al. (1993); Fochesato and Dias (2006)) intend to improve the efficiency of the BEM model. Finite difference methods
(see Engsig-Karup et al. (2009)) 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 (Wu and Taylor (1994); Ma and Yan (2006)) or a recent fully dispersive coupled-mode model described by Belibassakis et al. (2014). 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 Bingham et al. (2004) and Madsen et al. (2006) 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 Fenton and Rienecker (1982), 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 $\sigma$-transform model introduced by Chern et al. (1999) 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 Fructus et al. (2005) 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 Craig and Sulem (2001) for a flat bottom. This method was extended to a variable bottom by Guyenne and Nicholls (2005) and Craig et al. (2005) by introducing another operator depending only on the variation of the bottom. It was then improved numerically by Guyenne and Nicholls (2007).

In the present paper, we use the High-Order Spectral (HOS) method. This
non-linear potential method has been initially developed by West et al. (1987) and Dommermuth and Yue (1987) for a flat bottom, and extensively validated for different configurations in the LHEEA Lab., from regular waves up to irregular multidirectional wavefields (see Ducrozet et al. (2007, 2012b)). This model, named HOS-ocean is available as an open-source version<fnref n="1" />. The HOS formalism presents expansions identical to the DNO method, as demonstrated in Schäffer (2008), 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. Liu and Yue (1998) 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 Guyenne and Nicholls (2007), 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

<footnote n="1">¹https://github.com/LHEEA/HOS-ocean/wiki</footnote>
the original method presented hereafter has already been partially introduced
in Gouin et al. (2015). For documenting the accuracy of both schemes, a test
case with a constant variation of the bottom is computed to check the conver-
gence 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 nu-
merically 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.

## 2. Methods and algorithms

### 2.1. *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*₀ the mean depth and *β*(*x*) the bottom variation, such as −*h*(*x*) = −*h*₀ + *β*(*x*).
Thus, the domain considered is: −*h*₀ + *β*(*x*) ≤ *z* ≤ *η*(*x*) (see Fig.1).
<figure><img src="image_2.png" /><figcaption>Figure 1: Description of the fluid domain</figcaption></figure>

An infinite domain assumption is adopted with the assumption of periodic boundary conditions in the horizontal plane.

A potential flow formalism is used (incompressible and inviscid fluid, irrotational flow). Given these assumptions, the velocity $\vec{V}$ derives from a potential $\phi$: $\vec{V}(x, z, t) = \vec{\nabla}\phi$ and the continuity equation becomes the Laplace equation in the fluid domain D:

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

Following Zakharov (1968), both kinematic and dynamic non-linear free-surface boundary conditions (FSBC) are written in terms of surface quantities $\eta$ and $\tilde{\phi}$ ($\tilde{\phi}(x, t) = \phi(x, z = \eta, t)$), and expressed at the exact free-surface position $z = \eta$:

$$ \frac{\partial \eta}{\partial t} = \left( 1 + \left| \frac{\partial \eta}{\partial x} \right|^2 \right) \frac{\partial \phi}{\partial z} - \frac{\partial \tilde{\phi}}{\partial x} \cdot \frac{\partial \eta}{\partial x} \qquad (2) $$

$$ \frac{\partial \tilde{\phi}}{\partial t} = -g\eta - \frac{1}{2} \left| \frac{\partial \tilde{\phi}}{\partial x} \right|^2 + \frac{1}{2} \left( 1 + \left| \frac{\partial \eta}{\partial x} \right|^2 \right) \left( \frac{\partial \phi}{\partial z} \right)^2 \qquad (3) $$

To account for the time evolution of the quantities of interest $\eta$ and $\tilde{\phi}$ one only needs to evaluate the vertical velocity at the free surface $W(x, t) = \frac{\partial \phi}{\partial z}(x, z = \eta(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 :

$$ \frac{\partial \phi}{\partial x} \frac{\partial \beta}{\partial x} - \frac{\partial \phi}{\partial z} = 0 \quad \text{on} \quad z = -h_0 + \beta(x) \qquad (4) $$

To account for a varying bathymetry, Liu and Yue (1998) introduces an additional potential. Thus, the total potential $\phi_{tot}$ solution of the problem (which was named $\phi$ before) is expressed as:

$$ \phi_{tot} = \phi_{h_0} + \phi_{\beta} \qquad (5) $$

$\phi_{h_0}$ satisfies a Neumann condition on $z = -h_0$, therefore $\phi_{h_0}$ is solution of the problem at constant depth $h_0$:

$$ \frac{\partial \phi_{h_0}}{\partial z} (x, z = -h_0, t) = 0 \quad \text{on} \quad z = -h_0. \qquad (6) $$

$\phi_{\beta}$ allows the definition of the correct bottom boundary condition (Eq.4) and satisfies a Dirichlet condition on $z = 0$:

$$ \phi_{\beta} (x, z = 0, t) = 0 \quad \text{on} \quad z = 0. \qquad (7) $$

In 2D, the potentials are expanded on basis functions taking into account previous boundary conditions and the periodicity of the domain:

$$ \phi_{h_0} (x, z, t) = \sum_j A_j (t) \frac{\cosh(k_j (z + h_0))}{\cosh(k_j h_0)} e^{ik_j x} \qquad (8) $$

$$ \phi_{\beta} (x, z, t) = \sum_j B_j (t) \frac{\sinh(k_j z)}{\cosh(k_j h_0)} e^{ik_j x} \qquad (9) $$

with $k_j = j \frac{2\pi}{L_x}$ and $A_j$ and $B_j$ the modal amplitudes of $\phi_{h_0}$ and $\phi_{\beta}$ respectively.

## 2.2. High-Order-Spectral Method

The HOS model is a pseudo-spectral model initially developed in West et al. (1987); Dommermuth and Yue (1987). The potential is expressed as a truncated
power series of components $\phi_{tot}^{(m)}$ for $m = 1$ to $M$ ($M$ is the order of non-linearity 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).

$$ \tilde{\phi}_{tot} = \phi_{tot}(x, z = \eta, t) = \sum_{m=1}^{M} \sum_{n=0}^{M-m} \frac{\eta^n}{n!} \frac{\partial^n \phi_{tot}^{(m)}}{\partial z^n} (x, z = 0, t) \quad (10) $$

$$ \frac{\partial \phi_{tot}}{\partial z}(x, z = \eta, t) = W(x, t) = \sum_{m=1}^{M} \sum_{k=0}^{m-1} \frac{\eta^k}{k!} \frac{\partial^{k+1} \phi_{tot}^{(m-k)}}{\partial z^{k+1}}(x, 0, t) \quad (11) $$

Equations 11 can be reduced to Eq. (5.7) cited in Liu and Yue (1998) by an appropriate change of variables.

One more equation (the bottom condition) is needed to find the modal amplitudes of order $m: A_j^{(m)}(t)$ and $B_j^{(m)}(t)$. These modal amplitudes are the Fourier coefficients of the potentials Eqs. (8,9) at order $m$. Once $A_j^{(m)}(t)$ and $B_j^{(m)}(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 $\eta$ and $\tilde{\phi}_{tot}$ (as presented in Fig.3), which are used in a time-marching classical fourth-order Runge-Kutta scheme with an adaptive time step (see Fructus et al. (2005) 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 Ducrozet et al. (2007)). Note that a careful dealiasing is performed as explained in Bonnefoy et al. (2010) 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 Dommermuth and Yue (1987). The two following sections will present the addition due to the consideration of a non-uniform bottom.
<figure><img src="image_3.png" /><figcaption>Figure 2: Use of the bottom condition in the temporal solution.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Temporal solution of the FSBC.</figcaption></figure>
### 2.3. Original method

The way of using the bottom condition to find the modal amplitudes $B_j^{(m)}(t)$ (Eq.(5.4) from Liu and Yue (1998)) is depicted in Fig.2 and explained hereafter. The bottom condition (Eq.4) truncated at the order M reads:

$$ \sum_{m=1}^{M} \frac{\partial \phi_{tot}^{(m)}}{\partial x} \frac{\partial \beta}{\partial x} - \sum_{m=1}^{M} \frac{\partial \phi_{tot}^{(m)}}{\partial z} = 0 \quad \text{on } z = -h_0 + \beta(x) \quad (12) $$

By assuming that $\frac{\beta}{h_0} \ll 1$ we can write a Taylor expansion with respect to the mean depth $z = -h_0$:

$$ \sum_{m=1}^{M} \left[ \left( \frac{\partial \phi_{tot}^{(m)}}{\partial x} \frac{\partial \beta}{\partial x} \right)_{-h_0} + \beta \frac{\partial}{\partial z} \left( \frac{\partial \phi_{tot}^{(m)}}{\partial x} \frac{\partial \beta}{\partial x} \right)_{-h_0} + \dots + \frac{\beta^M}{M!} \frac{\partial^M}{\partial z^M} \left( \frac{\partial \phi_{tot}^{(m)}}{\partial x} \frac{\partial \beta}{\partial x} \right)_{-h_0} \right] \\ = \sum_{m=1}^{M} \left[ \left( \frac{\partial \phi_{tot}^{(m)}}{\partial z} \right)_{-h_0} + \beta \left( \frac{\partial^2 \phi_{tot}^{(m)}}{\partial z^2} \right)_{-h_0} + \dots + \frac{\beta^M}{M!} \left( \frac{\partial^{M+1} \phi_{tot}^{(m)}}{\partial z^{M+1}} \right)_{-h_0} \right] \quad (13) $$

Since we follow the potential theory (Eq.1) we have in two dimensions: $\frac{\partial^2 \phi_{tot}}{\partial z^2} = -\frac{\partial^2 \phi_{tot}}{\partial x^2}$, and we obtain:

$$ \sum_{m=1}^{M} \left[ \left( \frac{\partial \phi_{tot}^{(m)}}{\partial x} \frac{\partial \beta}{\partial x} \right)_{-h_0} + \dots + \frac{\beta^M}{(M)!} \frac{\partial^M}{\partial z^M} \left( \frac{\partial \phi_{tot}^{(m)}}{\partial x} \frac{\partial \beta}{\partial x} \right)_{-h_0} \right] \\ = \sum_{m=1}^{M} \left[ \left( \frac{\partial \phi_{tot}^{(m)}}{\partial z} \right)_{-h_0} - \beta \left( \frac{\partial^2 \phi_{tot}^{(m)}}{\partial x^2} \right)_{-h_0} - \dots - \frac{\beta^M}{M!} \left[ \frac{\partial^{M-1}}{\partial z^{M-1}} \left( \frac{\partial^2 \phi_{tot}^{(m)}}{\partial x^2} \right)_{-h_0} \right] \right] \quad (14) $$

We also assume $O(\beta) \equiv O(\frac{\partial \beta}{\partial x}) \equiv O(\eta)$ and we only keep terms of order $\eta^{(m)}$:

$$ \sum_{l=1}^{m-1} \frac{\beta^{l-1}}{(l-1)!} \frac{\partial^{l-1}}{\partial z^{l-1}} \left( \frac{\partial \phi_{tot}^{(m-l)}}{\partial x} \frac{\partial \beta}{\partial x} \right)_{-h_0} = \left( \frac{\partial \phi_{tot}^{(m)}}{\partial z} \right)_{-h_0} - \sum_{l=1}^{m-1} \frac{\beta^l}{l!} \left[ \frac{\partial^{l-1}}{\partial z^{l-1}} \left( \frac{\partial^2 \phi_{tot}^{(m-l)}}{\partial x^2} \right)_{-h_0} \right] \quad (15) $$

Thus we find the equations presented in Liu and Yue (1998):

$$ m = 1 : \frac{\partial \phi_{tot}^{(1)}}{\partial z} (x, -h_0, t) = 0 $$
$$m = 2, ..M : \frac{\partial \phi_{tot}^{(m)}}{\partial z}(x, -h_0, t) = \sum_{l=1}^{m-1} \frac{\partial}{\partial x} \left[ \frac{\beta^l}{l!} \frac{\partial^{l-1}}{\partial z^{l-1}} \left( \frac{\partial \phi_{tot}}{\partial x} \right)^{(m-l)} \right]_{z=-h_0} \quad (16)$$

Perturbation expansions are truncated at order M. This development (Eq.(16) with Eqs.(8, 9)) allows the computation of the modal amplitudes $B_{j_1}^{(m)}(t)$ at each order as functions of the $A_{j_1}^{(m)}(t)$. Modal amplitudes are independent of the position x, thus FFTs can still be used, preserving the numerical efficiency of the original HOS scheme.

## 2.4. Improved method

As described in Schäffer (2008), 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 Craig et al. (2005) and Guyenne and Nicholls (2007). Starting from Eq.(12) and assuming that $\frac{\beta}{h_0} \ll 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$:

$$\sum_{m=1}^{M} \left[ \left( \left( \frac{\partial \phi_{tot}^{(m)}}{\partial x} \frac{\partial \beta}{\partial x} \right)_{-h_0} + \beta \frac{\partial}{\partial z} \left( \frac{\partial \phi_{tot}^{(m)}}{\partial x} \frac{\partial \beta}{\partial x} \right)_{-h_0} + \dots + \frac{\beta^{M_b-1}}{(M_b-1)!} \frac{\partial^{M_b-1}}{\partial z^{M_b-1}} \left( \frac{\partial \phi_{tot}^{(m)}}{\partial x} \frac{\partial \beta}{\partial x} \right)_{-h_0} \right] \\ = \sum_{m=1}^{M} \left[ \left( \frac{\partial \phi_{tot}^{(m)}}{\partial z} \right)_{-h_0} - \beta \left( \frac{\partial^2 \phi_{tot}^{(m)}}{\partial x^2} \right)_{-h_0} - \dots - \frac{\beta^{M_b}}{(M_b)!} \left[ \frac{\partial^{M_b-1}}{\partial z^{M_b-1}} \left( \frac{\partial^2 \phi_{tot}^{(m)}}{\partial x^2} \right)_{-h_0} \right] \right] \quad (17)$$

In the left-hand side of Eq.17, the expansion is made up to order $M_b - 1$ because we assume that $O(\beta) \equiv O(\frac{\partial \beta}{\partial 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:

$$\sum_{m=1}^{M} \left( \frac{\partial \phi_{tot}^{(m)}}{\partial z} \right)_{-h_0} = \sum_{m'=1}^{M} \sum_{l=1}^{M_b} \left\{ \frac{\partial}{\partial x} \left[ \frac{\beta^l}{l!} \frac{\partial^{l-1}}{\partial z^{l-1}} \left( \frac{\partial \phi_{tot}^{(m')}}{\partial x} \right) \right]_{-h_0} \right\} \quad (18)$$

We adopt the hypothesis that each potential $\phi_{\beta}^{(m)}$ can be expressed as a truncated power series of components $\phi_{\beta}^{(m,l)}$ for $l = 0$ to $M_b$ (with $M_b$ the order of
non-linearity of the bottom):

$$ \phi_{\beta}^{(m)} = \sum_{l=1}^{M_b} \phi_{\beta}^{(m,l)} \qquad (19) $$

and the corresponding modal amplitudes are given by:

$$ B_j^{(m)} = \sum_{l=1}^{M_b} B_j^{(m,l)} \qquad (20) $$

Thus, at each order $m$, introducing $\phi_{tot}^{(m)} = \phi_{h_0}^{(m)} + \phi_{\beta}^{(m)}$ and keeping only the terms of maximum order $m + M_b$ (assuming as for the original method that $O(\beta) \equiv O(\frac{\partial \beta}{\partial x}) \equiv O(\eta)$) we get:

$$ \left( \frac{\partial \phi_{\beta}^{(m)}}{\partial z} \right)_{-h_0} = \sum_{l=1}^{M_b} \left( \frac{\partial \phi_{\beta}^{(m,l)}}{\partial z} \right)_{-h_0} = \sum_{l=1}^{M_b} \left\{ \left( \frac{\partial}{\partial x} \frac{\beta^l}{l!} \right) \left[ \frac{\partial^{l-1}}{\partial z^{l-1}} \left( \frac{\partial \phi_{h_0}^{(m)}}{\partial x} \right) \right]_{-h_0} \right\} + \sum_{l=1}^{M_b} \sum_{p=1}^{l-1} \left\{ \left( \frac{\partial}{\partial x} \frac{\beta^p}{p!} \right) \left[ \frac{\partial^{p-1}}{\partial z^{p-1}} \left( \frac{\partial \phi_{\beta}^{(m,l-p)}}{\partial x} \right) \right]_{-h_0} \right\} \quad (21) $$

We can then derive an iterative formula to calculate each $B_j^{(m,l)}$ at a given order $m$. This expression can be found in the DNO formalism in Guyenne and Nicholls (2007). Here we use the same formalism to compare our formulas with the ones given in their paper.

For $l$ odd:

$$ \begin{aligned} k B^{(m,l)} = & -k \left[ \frac{\beta^l}{l!} \operatorname{sech}(kh_0) k^{l-1} k \right] \\ & -k \left[ \sum_{p=2, \text{even}}^{l-1} \frac{\beta^p}{p!} k^{p-2} k k B^{(m,l-p)} - \sum_{p=1, \text{odd}}^{l-2} \frac{\beta^p}{p!} \tanh(kh_0) k^{p-2} k k F_{l-p} \right] \\ & +ik \frac{\partial \beta}{\partial x} \left[ \frac{\beta^{l-1}}{(l-1)!} \operatorname{sech}(kh_0) k^{l-2} k \right] \\ & +ik \frac{\partial \beta}{\partial x} \left[ \sum_{p=2, \text{even}}^{l-1} \frac{\beta^{p-1}}{(p-1)!} k^{p-3} k k B^{(m,l-p)} - \sum_{p=1, \text{odd}}^{l-2} \frac{\beta^{p-1}}{(p-1)!} \tanh(kh_0) k^{p-3} k k B^{(m,l-p)} \right] \end{aligned} \qquad (22) $$
And for $l$ even:

$$k B^{(m,l)} = -k \left[ \sum_{p=2,even}^{l-2} \frac{\beta^p}{p!} k^{p-2} k k B^{(m,l-p)} - \sum_{p=1,odd}^{l-1} \frac{\beta^p}{p!} \tanh(kh_0) k^{p-2} k k B^{(m,l-p)} \right] \\ + ik \frac{\partial \beta}{\partial x} \left[ \sum_{p=2,even}^{l-2} \frac{\beta^{p-1}}{(p-1)!} k^{p-3} k k B^{(m,l-p)} - \sum_{p=1,odd}^{l-1} \frac{\beta^{p-1}}{(p-1)!} \tanh(kh_0) k^{p-3} k k B^{(m,l-p)} \right] \quad (23)$$

If we compare Eq.(22,23) to the developments obtained by Guyenne and Nicholls (2007), one can notice that we obtain an additional part due to the gradient of the bottom (this part is missing in Guyenne and Nicholls (2007)). 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_j^{(m)}(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_j^{(m)}(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, \dots, M$. Section 3.1.3 shows the differences on the error made on the vertical velocity with the two methods.

## 2.5. Lateral boundary conditions and relaxation zones

With the present scheme periodic boundary conditions are required in the $x$-direction. 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 \le C_r(x) \le 1$ and an exact imposed solution $(\eta_{imp}, \tilde{\phi}_{imp})$ as described in Bingham and Zhang (2007) and Guyenne and Nicholls (2007). In our cases, the imposed solution will be either the stream solution of Rienecker and Fenton (1981) (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:

$$ \eta(x,t) = (1 - C_r(x)).\eta_{imp} + C_r(x).\eta(x,t) \quad (24) $$

$$ \tilde{\phi}(x,t) = (1 - C_r(x)).\tilde{\phi}_{imp} + C_r(x).\tilde{\phi}(x,t) \quad (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.

### 3. Regular wave propagation

A non-linear regular wave propagating over a flat bottom is considered. The initialisation of $\eta$ and $\tilde{\phi}$ is made by a non-linear regular wave solution of Rienecker and Fenton (1981) built at water depth $h$. The potential $\phi_{h0}$ is defined at the mean depth $h_0$ (which differs from $h$) and the potential $\phi_{\beta}$ must modify and correct the water depth $\beta$ ($\beta$ is a constant) which is introduced
artificially (cf Fig.4). The goal is to rebuild the potential $\phi_{tot}$ at the total water
depth $-h(x) = -h_0 + \beta(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.

### 3.1. Error on the vertical velocity.

In order to advance in time the unknowns on the free surface $\eta$ and $\tilde{\phi}$
(cf Eqs 2 and 3), we first need to evaluate the vertical velocity $W(x,t) =$
$\frac{\partial\phi}{\partial z}(x,z=\eta(x,t),t)$ associated to $\phi_{tot}$. This test case proposes to study the
error $\epsilon_w = \frac{W-W_{RF}}{W_{RF}}$ made on the evaluation of the vertical velocity once $\eta$ and
$\tilde{\phi}$ are known ($W_{RF}$ denotes the reference vertical velocity from Rienecker and
Fenton (1981)). The level of error correctly quantifies the accuracy of the model
in the propagation of such waves, as demonstrated in Bonnefoy et al. (2010).

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 $\frac{\beta}{h_0}$, 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 ($\frac{\beta}{h_0} = 0\%$) which was extensively validated
in Bonnefoy et al. (2010). 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 Gouin et al. (2014) 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.

#### 3.1.1. 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
<figure><img src="image_5.png" /><figcaption>Figure 4: Schematic view of the validation case</figcaption></figure>

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 Williams (1981) and the parametrized formula of Fenton (1990). 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 $\frac{\beta}{h_0}$. The lower surface is for $\frac{\beta}{h_0} = 0\%$ (i.e the initial HOS method) while the upper surface represents the convergence for a variation of the bottom $\frac{\beta}{h_0} = 25\%$. We observe that when $\frac{\beta}{h_0}$ 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 $M$-direction, 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 $\frac{\beta}{h_0} = 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
<figure><img src="image_6.png" /><figcaption>Figure 5: Error on the vertical velocity in function of N and M - kh=0.916 and ka=0.1.
Lower-surface : $\frac{\beta}{h_0}$ = 0%. Higher-surface : $\frac{\beta}{h_0}$ = 25%.</figcaption></figure>

more difficult to achieve in the *M*-direction). We also computed different combinations of the parameters *kh*, *ka*, $\frac{\beta}{h_0}$. An example of errors achieved with *N* = 64, *M* = 8 and different values of *kh* and *ka* is shown in Tab.1.

265

<figure><img src="image_7.png" /><figcaption>Figure 6: Error on the vertical velocity in function of N and M - kh=0.916 and ka=0.1.
Lower-surface : $\frac{\beta}{h_0}$ = 0% Higher-surface : $\frac{\beta}{h_0}$ = 50%.</figcaption></figure>

We conclude that with the original method, when *kh* increases (with *ka*
<table><caption>Table 1: Relative error on W, N = 64, M = 8.</caption><thead><tr><th></th><th colspan="3">kh=0.916</th></tr></thead><tbody><tr><th scope="row">β/h₀</th><td>k.a=0.01</td><td>k.a=0.05</td><td>k.a=0.1</td></tr><tr><th scope="row">0%</th><td>1.7 10<sup>-14</sup></td><td>4.0 10<sup>-10</sup></td><td>2.5 10<sup>-7</sup></td></tr><tr><th scope="row">5%</th><td>4.8 10<sup>-12</sup></td><td>1.8 10<sup>-9</sup></td><td>3.0 10<sup>-7</sup></td></tr><tr><th scope="row">25%</th><td>5.5 10<sup>-8</sup></td><td>2.6 10<sup>-6</sup></td><td>3.4 10<sup>-5</sup></td></tr><tr><th scope="row">50%</th><td>3.9 10<sup>-6</sup></td><td>1.1 10<sup>-4</sup></td><td>1.2 10<sup>-3</sup></td></tr></tbody><thead><tr><th></th><th colspan="3">ka=0.1</th></tr></thead><tbody><tr><th scope="row">β/h₀</th><td>k.h=0.67</td><td>k.h=3</td><td>k.h=10</td></tr><tr><th scope="row">0%</th><td>5.7 10<sup>-8</sup></td><td>6.3 10<sup>-9</sup></td><td>6.2 10<sup>-9</sup></td></tr><tr><th scope="row">5%</th><td>2.2 10<sup>-7</sup></td><td>6.3 10<sup>-9</sup></td><td>6.2 10<sup>-9</sup></td></tr><tr><th scope="row">25%</th><td>3.9 10<sup>-5</sup></td><td>1.8 10<sup>-6</sup></td><td>2.6 10<sup>-7</sup></td></tr><tr><th scope="row">50%</th><td>1.1 10<sup>-3</sup></td><td>3.8 10<sup>-4</sup></td><td>1.8 10<sup>-4</sup></td></tr></tbody></table>

and $\frac{\beta}{h_0}$ fixed) the convergence rate is better. This result was expected because when $\frac{h}{\lambda} > 2$ (i.e. $kh > \pi$), 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 $\frac{\beta}{h_0}$ increase (with the other parameters fixed), the convergence rate decreases. Of course there is a limit in the admissible values of $ka$ and $\frac{\beta}{h_0}$. The steepness can not exceed the theoretical highest stable wave given by the formula of Fenton (1990) and $\frac{\beta}{h_0}$ 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.

### 3.1.2. 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 $\frac{\beta}{h_0}$.

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
$\frac{\beta}{h_0} = 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 be drawn with increasing $kh$, $ka$ and $\frac{\beta}{h_0}$ and comparisons between the two approaches are provided in the following section.

290

<figure><img src="image_8.png" /><figcaption>Figure 7: Error on the vertical velocity in function of $M$ and $M_b - kh = 0.916$, $ka = 0.1$ and $\frac{\beta}{h_0} = 25\%$.</figcaption></figure>

### 3.1.3. Comparison of the two methods

295

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 $\frac{\beta}{h_0}$.

300

Fig.8 and Fig.9 present the error achieved on the vertical velocity as function of $M$ and $M_b$ for $\frac{\beta}{h_0} = 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
(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.

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_j^{(m)}(t)$ 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.

The influence of $\frac{\beta}{h_0}$ is represented on Fig.10 and Fig.11. In Fig.10, for $M_b \ge 3$ the error achieved on the vertical velocity is converged so the corresponding curves are superimposed. For these cases we see that for a small variation of the bottom (Fig.10), both methods present the same accuracy. But when $\frac{\beta}{h_0}$ increases, the improved method seems much more accurate, for the same reason as before: when $\frac{\beta}{h_0}$ 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 $\frac{\beta}{h_0} = 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 \ge 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 Rienecker and Fenton (1981). 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 Bonnefoy et al. (2010)). It happens even if the bottom is flat ($ie \frac{\beta}{h_0} = 0\%$), and we have checked in Gouin et al. (2015) that neither method induces a prohibitive phase shift.

<figure><img src="image_9.png" /><figcaption>Figure 8: Error on the vertical velocity as function of $M$ and $M_b - \frac{\beta}{h_0} = 40\% - kh = 0.916$ and $ka = 0.01$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Error on the vertical velocity as function of $M$ and $M_b - \frac{\beta}{h_0} = 40\% - kh = 0.916$ and $ka = 0.15$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Error on the vertical velocity as function of M and $M_b - \frac{\beta}{h_0} = 10\% - kh = 0.916$ and $ka = 0.15$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Error on the vertical velocity as function of M and $M_b - \frac{\beta}{h_0} = 75\% - kh = 0.916$ and $ka = 0.15$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Error on the vertical velocity as function of M and $M_b - \frac{\beta}{h_0} = 40\% - kh = 0.5$ and $ka = 0.1$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Error on the vertical velocity as function of M and $M_b - \frac{\beta}{h_0} = 40\% - kh = 3$ and $ka = 0.1$.</figcaption></figure>
## 3.2. 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^{\frac{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) \qquad (26)$$

While for the improved method this number is:

$$n = 4M_b M \qquad (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 < \frac{M-1}{2} \quad \text{or} \quad M < \frac{M_b-1}{2} \qquad (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^{\frac{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 $\frac{\beta}{h_0} = 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 $\frac{\beta}{h_0} = 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.

<figure><img src="image_15.png" /><figcaption>Figure 14: Computational effort with respect to $N$. $M = 8$.</figcaption></figure>

### 3.3. *Partial conclusion*

Both methods can compute various cases depending on the value of $\frac{\beta}{h_0}$, $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.
390

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.

395

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.

# 4. Validation

400

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.

405

## 4.1. Bragg reflection

410

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 Davies and Heathershaw (1984) and numerically by Bingham and Zhang (2007); Guyenne and Nicholls (2007); Liu and Yue (1998).

415

For small incident waves and small bottom slopes, reflection near Bragg resonances is well predicted by the multiple-scale perturbation theory of Mei
(1985). Here we analyse non-linear effects using the original method. Note that the improved method gives exactly the same results. The conditions of Davies and Heathershaw (1984) are used to compare with experiments.

420

<figure><img src="image_16.png" /><figcaption>Figure 15: Bottom topography with a patch of 10 sinusoidal ripples of amplitude $d = 0.1$ and slope $k_b d = 0.31$.</figcaption></figure>

The bottom ripple patch is defined as the variation around the mean water depth $h = h_0 + \beta(x)$ by :

$$ \beta(x) = d \sin(k_b x) \quad \text{for} \quad \frac{-L_0}{2} \le x \le \frac{L_0}{2} \qquad (29) $$

as depicted in Fig.15. The free surface is located at $z = 0$m and the mean water depth is at $h_0 = 1$m. The ripple slope is $k_b d = 0.31$, the ripple amplitude is $d = 0.1$m and the length of the patch is $\frac{L_0}{\lambda_b} = 10$ (i.e. we consider a patch of 10 sinusoidal ripples of wavelength $\lambda_b = \frac{2\pi}{k_b}$). The incident wave is set at the linear resonance condition of $k = \frac{k_b}{2}$ with a wave steepness $ka = 0.05$.

425

We perform simulations with $\frac{N}{\lambda} = 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 Liu and Yue (1998). The local reflection coefficient $R(x)$ is then extracted

430
using the method of Suh et al. (2001) and analysed hereafter.

The results are presented in Fig.16 along with the experimental measure-
ments of Davies and Heathershaw (1984) and the solution given by the pertur-
bation theory of Mei (1985). Ripples are located in $x/\lambda_b \in [-5; 5]$. It appears

<figure><img src="image_17.png" /><figcaption>Figure 16: Bragg reflection from a sinusoidal bottom ripple patch over $-5\lambda_b \le x \le 5\lambda_b$ with $ka = 0.05$ and $k_b d = 0.31$.</figcaption></figure>

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 Liu and Yue
(1998). Thus the proposed methods are validated for a classical configuration
with a small variation of the bottom ($\frac{\beta}{h_0} = 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.

4.2. 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 Bingham and Agnon (2005); Guyenne and Nicholls (2007).

450

We used a bottom profile defined by:

$$ \beta(x) = \frac{h_0 - h_1}{2} \left[ 1 + \tanh \left( \frac{\sin(\pi x/L)}{1 - (2x/L)^2} \right) \right], \quad -\frac{L}{2} \le x \le \frac{L}{2} \quad (30) $$

<figure><img src="image_18.png" /><figcaption>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 $\lambda = 2$ m over a shoaling bottom profile at steady-state.</figcaption></figure>

455

with $h_0 = 0.5$ m, $h_1 = 0.1$ m and $L = 24$ m. Here $\frac{\beta}{h_0} = 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 $\lambda_0 = 2$ m 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).

460
wavelengths) and a second absorbing zone (10 wavelengths) (see Fig.17).

465

The mean water depth $h_0$ is chosen at 0.3 m to reduce the maximum $\frac{\beta}{h_0}$ (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:

$$ \frac{a}{a_0} = \left[ \frac{k \left( 1 + 2h_0 k_0 / \sinh(2h_0 k_0) \right)}{k_0 \left( 1 + 2hk / \sinh(2hk) \right)} \right]^{1/2} \quad (31) $$

<figure><img src="image_19.png" /><figcaption>Figure 18: Comparison between the envelope predicted by the linear theory and the envelope obtained with the computed wave profile.</figcaption></figure>

470

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 Bingham and Agnon (2005)). Thus this test case validates the improved method even for non-negligible relative bottom variations $\frac{\beta}{h_0}$.

475
## 5. 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.

### 5.1. Mild slope

This test case has been studied both experimentally by Dingemans (1994); Luth et al. (1994) and numerically by Guyenne and Nicholls (2007); Bingham and Zhang (2007). The bottom variation is defined by:

$$ \beta(x) = \begin{cases} 0.05(x-6) & \text{for } 6 \le x \le 12, \\ 0.3 & \text{for } 12 \le x \le 14, \\ 0.3 - 0.1(x-14) & \text{for } 14 \le x \le 17, \\ 0 & \text{elsewhere,} \end{cases} \qquad (32) $$

and can be seen in Fig.19. It has been scaled with a factor of two in comparison with the experiment of Dingemans (1994). 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 $\frac{\beta}{h_0} = 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 Rienecker and Fenton (1981) 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.
<figure><img src="image_20.png" /><figcaption>Figure 19: Sketch of the submerged bar as described in the experiments of Dingemans (1994). Blue markers correspond to the gauges in the experiments.</figcaption></figure>

With the two methods the convergence and steady-state are reached after $t = 25T$ with a number of nodes per wavelength $\frac{N}{\lambda} = 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.

### 5.1.1. 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 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><img src="image_21.png" /><figcaption>Figure 20: Snapshot of the surface elevation at a steady-state (scaled by a factor of 3).</figcaption></figure>

The comparison between our numerical results and the experimental data is very good, and similar to the results obtained with other numerical methods (Guyenne and Nicholls (2007); Bingham and Zhang (2007)). 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.

### 5.1.2. Harmonic analysis

For a deeper comparison, a harmonic analysis of the surface elevation is run as presented in Bingham and Zhang (2007). 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 Bingham and Zhang (2007), 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 $\frac{\beta}{h_0} = 0.75$ over
<figure><img src="image_22.png" /><figcaption>Figure 21: Time series of measured and computed surface elevations at 6 positions. Solid line: calculations. Points: measurements Luth et al. (1994).</figcaption></figure>
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).

<figure><img src="image_23.png" /><figcaption>Figure 22: Harmonic analysis at steady-state for the mild slope case. $\frac{N}{\lambda} = 40$, $M = 3$ and $M_b = 16$.</figcaption></figure>

5.2. 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 Ohyama et al. (1995).
In this case, $\frac{\beta}{h_0} = 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 Shen
et al. (2004); Gobbi and Kirby (1998); Benxia and Xiping (2009) as a validation
case for their numerical methods.
The bottom variation is presented in Fig.23 and defined for this steeper case by:

$$ \beta(x) = \begin{cases} 0.5(x - 2.6) & \text{for } 2.6 \le x \le 3.3, \\ 0.35 & \text{for } 3.3 \le x \le 4.8, \\ 0.35 - 0.5(x - 4.8) & \text{for } 4.8 \le x \le 5.5, \\ 0 & \text{elsewhere,} \end{cases} \quad (33) $$

<figure><img src="image_24.png" /><figcaption>Figure 23: Sketch of the submerged bar from Ohyama et al. (1995).</figcaption></figure>

The initialisation of $\eta$ and $\tilde{\phi}$ is based on a solution of Rienecker and Fenton (1981) 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/\lambda = 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
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 Ohyama et al. (1995).

565

<figure><img src="image_25.png" /><figcaption>Figure 24: Time series of measured and computed surface elevations at 2 positions. Solid line: calculations. Points: measurements Ohyama et al. (1995).</figcaption></figure>

The comparison between our numerical results and the experimental data
is also good even if we observe slight discrepancies with respect to the exper-
iments. In this steeper case, the maximum variation of the bottom is almost
the same as in the mild slope case (β/h₀ = 70%) but the gradient is much larger.
Nevertheless, the large slopes are more difficult to simulate with the proposed
models. However, in this extreme configuration, we still manage to obtain rea-
sonable 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

570
presented in Shen et al. (2004) and almost the same as in Gobbi and Kirby (1998); Benxia and Xiping (2009).

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 these 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.

## 6. 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 Liu and Yue (1998), while the second method is an extension of the work of Guyenne and Nicholls (2007) 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_{dlog}(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
(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 Liu and Yue (1998); Bingham and Zhang (2007); Guyenne and Nicholls (2007). 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 Bingham and Agnon (2005); Guyenne and Nicholls (2007).

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
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.

References

Belibassakis, K., Athanassoulis, G., Gerostathis, T., 2014. Directional wave spectrum transformation in the presence of strong depth and current inhomogeneities by means of coupled-mode model. Ocean Engineering 87, 84 – 96.

Benxia, L., Xiping, Y., 2009. Wave decomposition phenomenon and spectrum evolution over submerged bars. Acta Oceanologica Sinica 29.

Bingham, H., Fuhrman, D., Jamois, E., Kimmoun, O., 2004. Nonlinear wave interaction with bottom-mounted structures by a high-order Boussinesq method, in: 19th International Workshop on Water Waves and Floating Bodies, Proceedings, pp. 9–12.

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

Bingham, H.B., Agnon, Y., 2005. A Fourier-Boussinesq method for nonlinear water waves. European Journal of Mechanics - B/Fluids 24, 255–74.

Bonnefoy, F., Ducrozet, G., Le Touze, D., Ferrant, P., 2010. Time-domain simulation of nonlinear water waves using spectral methods. Advances in Numerical Simulation of Nonlinear Water Waves 11.

Chern, M., Borthwick, A., Taylor, R.E., 1999. A pseudospectral σ-transformation model of 2-D nonlinear waves. Journal of Fluids and Structures 13, 607–30.
Craig, W., Guyenne, P., Nicholls, D.P., Sulem, C., 2005. Hamiltonian long-wave expansions for water waves over a rough bottom, The Royal Society. pp. 839-73.

Craig, W., Sulem, C., 2001. Numerical simulation of gravity waves. Journal of Computational Physics 108, 73-83.

Davies, A., Heathershaw, A., 1984. Surface-wave propagation over sinusoidally varying topography. Journal of Fluid Mechanics 144, 419-43.

Dingemans, M., 1994. Comparison of computations with Boussinesq-like models and laboratory measurements. Delft Hydraulics memo H1684.12.

Dommermuth, D., Yue, D., 1987. A high-order spectral method for the study of non-linear gravity waves. Journal of Fluid Mechanics 184, 267-88.

Ducrozet, G., Bingham, H., Engsig-Karup, A.P., Bonnefoy, F., Ferrant, P., 2012a. A comparative study of two fast nonlinear free-surface water wave models. International Journal for Numerical Methods in Fluids 69, 1818-34.

Ducrozet, G., Bonnefoy, F., Le Touze, D., Ferrant, P., 2007. 3-D HOS simulations of extreme waves in open seas. Natural Hazards and Earth System Sciences 7, 109-22.

Ducrozet, G., Bonnefoy, F., Le Touze, D., Ferrant, P., 2012b. A modified high-order spectral method for wavemaker modeling in a numerical wave tank. European Journal of Mechanics - B/Fluids 34, 19-34.

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

Fenton, J., 1990. Nonlinear wave theories. The Sea, Vol.9: Ocean Engineering Science. pp. 3-25.
Fenton, J., Rienecker, M., 1982. A Fourier method for solving nonlinear water-wave problems: application to solitary-wave interactions. Journal of Fluid Mechanics 118, 411–43.

Fochesato, C., Dias, F., 2006. A fast method for nonlinear three-dimensional free-surface waves. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 462, 2715–35.

Fructus, D., Clamond, D., Grue, J., Kristiansen, O., 2005. An efficient model for three-dimensional surface wave simulations. Journal of Computational Physics 205, 665–85.

Gobbi, M.F., Kirby, J.T., 1998. Wave evolution over submerged sills: Tests of a high-order Boussinesq model. Coastal Engineering 37, 57–96.

Gouin, M., Ducrozet, G., Ferrant, P., 2014. Development of a highly nonlinear model for wave propagation over a variable bathymetry, in: 29th international workshop on water waves and floating bodies, Proceedings.

Gouin, M., Ducrozet, G., Ferrant, P., 2015. Development and validation of a highly nonlinear model for wave propagation over a variable bathymetry, in: 34th International Conference on Ocean, Offshore and Arctic Engineering, Proceedings.

Grilli, S., Guyenne, P., Dias, F., 2001. A fully non-linear model for three-dimensional overturning waves over an arbitrary bottom. International Journal for Numerical Methods in Fluids 35, 829–67.

Guyenne, P., Nicholls, D., 2005. Numerical simulation of solitary waves on plane slopes. Mathematics and Computers in Simulation 69, 269–81.

Guyenne, P., Nicholls, D., 2007. A high-order spectral method for nonlinear water waves over moving bottom topography. Journal of Scientific Computing 30, 81–101.
Korsmeyer, F.T., Yue, D.K.P., Nabors, K., White, J., 1993. Multipole-accelerated preconditioned iterative methods for three-dimensional potential problems, in: 15th international Conference on Boundary Element, Proceedings.

Liu, Y., Yue, D., 1998. On generalized Bragg scattering of surface waves by bottom ripples. Journal of Fluid Mechanics 356, 297-326.

Luth, H., Klopman, G., Kitou, N., 1994. Kinematics of waves breaking partially on an offshore bar. Delft Hydraulics memo H1573.

Ma, Q., Yan, S., 2006. Quasi ALE finite element method for nonlinear water waves. Journal of Computational Physics 212, 52-72.

Madsen, P.A., Fuhrman, D.R., Wang, B., 2006. A Boussinesq-type method for fully nonlinear waves interacting with a rapidly varying bathymetry. Coastal Engineering 53, 487-504.

Mei, C.C., 1985. Resonant reflection of surface water waves by periodic sandbars. Journal of Fluid Mechanics 152, 315-35.

Ohyama, T., Kioka, W., Tada, A., 1995. Applicability of numerical models to nonlinear dispersive waves. Coastal Engineering 24, 297-313.

Rienecker, M., Fenton, J., 1981. A Fourier approximation method for steady water waves. Journal of Fluid Mechanics 104, 119-37.

Schäffer, H., 2008. Comparison of Dirichlet-Neumann operator expansions for nonlinear surface gravity waves. Coastal Engineering 55, 288-94.

Shen, Y., Ng, C., Zheng, Y., 2004. Simulation of wave propagation over a submerged bar using the VOF method with a two-equation $k-\epsilon$ turbulence modeling. Ocean Engineering 31, 87-95.

Suh, K.D., Park, W.S., Park, B.S., 2001. Separation of incident and reflected waves in wavecurrent flumes. Coastal Engineering 43, 149-59.
West, B., Brueckner, K., Janda, R.S., Milder, D.M., Milton, R.L., 1987. A new numerical method for the surface hydrodynamics. Journal of Geophysical Research 92, 11803-24.

740 Westphalen, J., Greaves, D., Raby, A., Hu, Z., Causon, D., Mingham, C., Omidvar, P., Stansby, P., D. Rogers, B., 2014. Investigation of wave-structure interaction using state of the art CFD techniques. Open Journal of Fluid Dynamics 4, 18-43.

745 Williams, J.M., 1981. Limiting gravity waves in water of finite depth. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 302, 139-88.

Wu, G.X., Taylor, R.E., 1994. Finite element analysis of two-dimensional non-linear transient water waves. Applied Ocean Research 16, 363-72.

Zakharov, V., 1968. Stability of periodic waves of finite amplitude on the surface of a deep fluid. Journal of Applied Mechanics and Technical Physics 9, 190-4.

## Vitae

Maïté Gouin is a PhD student at the LHEEA Lab. of Ecole Centrale Nantes (ECN) since February, 2013. The research is sponsored by the IRT Jules Verne (Research Institute of Technology) and deals with the development of a non-linear spectral model for wave propagation over a variable bathymetry.

Guillaume Ducrozet is assistant professor at the LHEEA Lab. of Ecole Centrale Nantes (ECN) since 2010. He is in charge of the developments of numerical models for non-linear gravity waves in oceans and wave tanks. He works on the development of efficient numerical models in this concern and especially with pseudo-spectral methods. He graduated from ECN with a Master of Engineering in free surface hydrodynamics in 2004, followed by a PhD obtained in 2007.
His PhD was dealing with the modeling of non-linear processes during genera-
tion and propagation of gravity waves.

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 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.

770