2.1. High-order spectral method

We consider a 3D rectangular fluid domain, periodic in the horizontal directions and equipped with a Cartesian coordinate system, as shown in figure 1. Here $z=\eta (x,y,t)$ denotes the free surface elevation and $z=0$ is the still water level. We use $z=-h_0+\beta (x,y)$ to denote the bottom topography, where $h_0$ is the mean water depth and $\beta (x,y)$ is the bottom variation.

Figure 1. (a) Smooth functions $\chi _{d}$ for damping zones (solid line) and $\chi _{nl}$ for spatial nonlinear adjustment (dashed line). (b) Computational domain with damping zones and influx line (red) for wave generation.

A potential theory is used assuming inviscid fluid and irrotational flow. The fluid velocity is expressed by the velocity potential $\phi$, $\boldsymbol {V}=\boldsymbol {\nabla }\phi$. Assuming the fluid is incompressible, the continuity equation becomes the Laplace equation

(2.1)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\nabla} \phi = 0. \end{equation}

We consider periodic boundary conditions in the horizontal plane. The bottom boundary condition is

(2.2)\begin{equation} \boldsymbol{\nabla}\phi \boldsymbol{\cdot} \boldsymbol{\nabla}\beta - \frac{\partial\phi}{\partial z} = 0 \quad \text{at}\ z={-}h_0+\beta(x,y). \end{equation}

As pointed out in Zakharov (Reference Zakharov1968), the kinematic and dynamic nonlinear free surface boundary conditions can be written in terms of surface elevation $\eta$ and surface potential $\phi _s=\phi (x,y,z=\eta,t)$

(2.3)\begin{equation} \left.\begin{gathered} \frac{\partial \eta}{\partial t} = W \left( 1+ \left\vert \boldsymbol{\nabla} \eta \right\vert^2\right) - \boldsymbol{\nabla}\phi_s \boldsymbol{\cdot} \boldsymbol{\nabla}\eta \\ \frac{\partial\phi_s}{\partial t} ={-}g \eta -\frac{1}{2}\left\vert\boldsymbol{\nabla}\phi_s\right\vert^2 + \frac{1}{2}W^2\left(1+\left\vert\boldsymbol{\nabla}\eta\right\vert^2\right), \end{gathered}\right\} \end{equation}

where $g$ is the acceleration of gravity and $W={\partial \phi }/{\partial z}\vert _{z=\eta }$ is the vertical velocity at the surface. The surface vertical velocity $W$ needs to be evaluated to solve the system. For flat bottom, an efficient pseudo-spectral method to calculate the surface vertical velocity $W$ in terms of $\eta$ and $\phi _s$ was initially introduced by Dommermuth & Yue (Reference Dommermuth and Yue1987) and West et al. (Reference West, Brueckner, Janda, Milder and Milton1987), and it has been known as the high-order spectral method (HOSM). Extension to variable depth was discussed in Gouin, Ducrozet & Ferrant (Reference Gouin, Ducrozet and Ferrant2016,  Reference Gouin, Ducrozet and Ferrant2017). A brief summary of the HOSM for variable depth by Gouin et al. (Reference Gouin, Ducrozet and Ferrant2017) is presented in the following.

The velocity potential $\phi$ is expressed as a truncated power series $\phi =\sum _{m=1}^M \phi ^{(m)}$ where $m$ is the nonlinear order of $\phi ^{(m)}$ and $M$ is the nonlinear order of the method which can be freely chosen. The velocity potential evaluated at the free surface $\phi _s$ is then expanded in a Taylor series around the still water level $z=0$.

For an uneven bottom, an additional velocity potential $\phi _B$ is expressed as a truncated power series $\phi _B=\sum _{m=1}^M \phi _B^{(m)}$ and $\phi _B^{(m)}=\sum _{l=1}^{M_B}\phi _B^{(m,l)}$ where $M_B$ is the order for the bottom boundary condition which can be different from $M$ and can be freely chosen. Thus, the velocity potential $\phi ^{(m)}$ is expressed as the sum of two components

(2.4)\begin{equation} \left.\begin{gathered} \phi^{(m)} = \phi_0^{(m)} + \phi_B^{(m)}\\ \phi_0^{(m)}(x,y,z,t) = \sum_{p,q} A_{pq}^{(m)}(t) \frac{\cosh(k_{pq}(z+h_0))}{\cosh(k_{pq}h_0)} \exp({{\rm i}(k_{x_p}x +k_{y_q}y)})\\ \phi_B^{(m)}(x,y,z,t) = \sum_{p,q} B_{pq}^{(m)}(t) \frac{\sinh(k_{pq}z)}{\cosh(k_{pq}h_0)} \exp({{\rm i}(k_{x_p}x +k_{y_q}y)}) \end{gathered}\right\} \end{equation}

with $k_{x_p}=(p-{N_x}/{2}+1)({2{\rm \pi} }/{L_x})$ for $p=0,\ldots,N_x-1$, $k_{y_q}=(q-{N_y}/{2}+1)({2{\rm \pi} }/{L_y})$ for $q=0,\ldots,N_y-1$ and $k_{pq}=\sqrt {k_{x_p}^2+k_{y_q}^2}$. Here $N_x$ and $N_y$ are the number of wave components in the $x$ and $y$ directions, respectively. We choose $N_x$ and $N_y$ to be even numbers. We use $A_{pq}^{(m)}(t)$ and $B_{pq}^{(m)}(t)=\sum _{l=1}^{M_B}B_{pq}^{(m,l)}(t)$ to denote the modal amplitudes of $\phi _0^{(m)}$ and $\phi _B^{(m)}$, respectively. With the velocity potential in (2.4) and Taylor expansion of the bottom boundary condition in (2.2), we can calculate each $B_{pq}^{(m,l)}(t)$.

For evaluation of the vertical velocity at the surface, $W$, a similar expansion series of surface potential is applied for the vertical velocity at the surface $W=\sum _{m=1}^{M}W^{(m)}$. Then, Taylor expansion around the still water level $z=0$ leads to a triangular system for $W^{(m)}$ to be solved iteratively. Once the vertical velocity at the surface $W^{(m)}$ has been computed, then the free surface elevation $\eta$ and surface potential $\phi _s$ can be marched in time through (2.3).

This method is based on a Taylor expansion of the bottom boundary condition with respect to the mean water depth and solved numerically by pseudo-spectral method. The bottom and its gradient need to be continuous to avoid any instabilities. For larger bottom variation $\beta /h_0$, the error on the vertical velocity at the surface, $W$, becomes larger and the convergence is slower but still exists by increasing number of wave components and orders of the method $M$ and $M_B$. A detailed convergence study in one horizontal dimension with respect to the number of wave component $N_x$ and orders $M$ and $M_B$ is presented for variable depth in Gouin et al. (Reference Gouin, Ducrozet and Ferrant2016).

2.2. Wave generation and damping zones

For a flat bottom, the simulation can be initiated by specifying the surface elevation $\eta$ and the surface potential $\phi _s$ at initial time $t=0$. The prescribed initial condition has to be periodic because of the pseudo-spectral method. For an uneven bottom, it is convenient to start simulations from rest ($\eta =\phi _s=0$) and use a wave generator inside the domain. Embedded wave generation is implemented with spatial nonlinear adjustment as described in Lie, Adytia & van Groesen (Reference Lie, Adytia and van Groesen2014). To ensure the periodicity, a couple of damping zones of length $\lambda _d$ are also applied.

The dynamic equations with embedded wave generation and damping zones become

(2.5)\begin{equation} \left.\begin{gathered} \frac{\partial \eta}{\partial t} = RHS_{lin} + RHS_{nonlin}\chi_{nl} - c_d\eta\chi_d + S_{x_0}\\ \frac{\partial \phi_s}{\partial t} = RHS_{lin} + RHS_{nonlin}\chi_{nl} - c_d\phi_s\chi_d \end{gathered}\right\} \end{equation}

where $RHS_{lin}$ and $RHS_{nonlin}$ are the linear and nonlinear parts of right-hand sides of the dynamic equations, $c_d$ is the damping coefficient, $S$ is the source term for embedded wave generation and $\chi _{nl}$ and $\chi _d$ are smooth functions. Illustration of the chosen functions for nonlinear adjustment and damping zones is shown in figure 1. In the damping zone, the solutions decay exponentially as $\textrm {e}^{-c_dT_d}$, where $T_d$ is the travel time of the wave in the damping zone. The value of the damping coefficient $c_d$ and the length of the damping zone $\lambda _d$ are coupled. As an example, for exponential decay of $10^{-3}\approx \textrm {e}^{-7}$, the length of the damping zone should be at least $\lambda _d=7V_g/c_d$, where $V_g$ is the group velocity.

2.3. Wave kinematics calculation

The HOSM calculates the surface vertical velocity $W$ in terms of surface elevation $\eta$ and surface potential $\phi _s$. For wave kinematics calculation, it is convenient to use (2.4) directly after the modal amplitudes $A_{pq}$ and $B_{pq}$ are known. We refer this method as the direct method in this paper. We consider a solitary wave and Stokes waves on a flat bottom to check the validity of the direct method.

Solitary waves with amplitude over depth $A/h=0.4$ and $A/h=0.6$ are chosen. We took the surface elevation and the surface potential for the solitary waves from the exact solution of the Euler equations by Dutykh & Clamond (Reference Dutykh and Clamond2014) as input to calculate the interior horizontal velocity with the direct method. Figure 2 shows the comparison of the horizontal velocity profile of solitary wave at the crest by the direct method and the reference solution from Dutykh & Clamond (Reference Dutykh and Clamond2014). For a rather extreme solitary wave $A/h=0.6$, the direct method gives inaccurate horizontal velocity close to the surface. For a milder case $A/h=0.4$, the accuracy of the direct method can be improved by increasing the nonlinear order.

Figure 2. Horizontal velocity profiles of solitary waves with (a) $A/h=0.4$ and (b) $A/h=0.6$.

We consider a Stokes waves in deep water ($kh=10$) with $ka=0.1$ and $ka=0.3$, where $a$ is the first-order amplitude, $k$ is the wavenumber and $h$ is the water depth. We use the fifth-order solution by Fenton (Reference Fenton1985) to obtain the surface elevation and surface potential. The surface potential is calculated from the fifth-order velocity potential at $z=\eta$, retaining fifth-order terms. Figure 3 shows the comparison of horizontal velocity under the wave crest by the direct method with Fenton's solution. For Stokes wave with steepness $ka=0.1$, the direct method gives a good agreement compared with Fenton's solution. Meanwhile, the horizontal velocity by the direct method starts to deviate from Fenton's solution close to the surface when the steepness is increased to $ka=0.3$.

Figure 3. Horizontal velocity profiles of Stokes waves with (a) $ka=0.1$ and (b) $ka=0.3$.

From both test cases, we see that the direct method is inaccurate to calculate the horizontal velocity near the surface for strongly nonlinear cases. However, the direct method is still valid down to the bottom and for low steepness parameter $A/h$ in shallow water or $ka$ in deep water. This problem has been pointed out in Bateman, Swan & Taylor (Reference Bateman, Swan and Taylor2003) where the $H$ and $H_2$ operators were introduced. For flat bottom, the direct method is equivalent to the $H$ operator method. The $H_2$ operator method was proposed in Bateman et al. (Reference Bateman, Swan and Taylor2003) as a remedy. For varying bottom, the more advanced variational Boussinesq model was proposed to calculate the wave kinematics in Lawrence, Gramstad & Trulsen (Reference Lawrence, Gramstad and Trulsen2021a).

In the present paper, we focus on the interior horizontal fluid velocity where it is not close to the water surface or the bottom. Therefore, we still use the direct method because it is efficient and sufficiently accurate to study the evolution of the statistical properties of velocity field.