Waves Generated by the Horizontal Motions of a Bottom Disturbance

Abstract

Waves generated by a horizontally moving disturbance on the seabed have been studied by developing two numerical models, namely, the Navier–Stokes and the Green–Naghdi equations. Various geometries of the bottom disturbances are considered, and waves generated due to a single motion and multiple oscillatory motions of the bottom disturbances are investigated by the two models. Discussion is provided on how the motion of the disturbance on the seafloor results in the generation of surface waves. The wave-field parameters investigated include the surface elevation, velocity, pressure fields and wave celerity. A parametric study is conducted to assess the effect of the geometry of the disturbance and the kinematic characteristics on the wave generation. It is shown that both linear and nonlinear waves can be generated by a horizontally moving disturbance on the seabed. Long waves, followed by a series of dispersive waves, are produced by the single motion of the bottom disturbance. It is also found that, under appropriate conditions, there would be a balance between nonlinearity and dispersion, such that the generated waves propagate over a flat seafloor with little to no change in their form and shape.

1. Introduction

Many geographical tsunamis are generated due to moving submarine landslides or earthquakes on the seafloor—see, e.g., [1,2,3]. Submarine landslides can assume arbitrary geometrical shapes with different movements—see, e.g., [4,5,6]. Recreating the wave generation by a moving disturbance on the seafloor to accurately recreate the generation and propagation of tsunami is of great interest.

Solitary wave and solitary-like waves are used commonly to model tsunamis—see, e.g., [7,8,9,10,11]. However, it has been found, e.g., by [12], that the solitary wave theory is not necessarily appropriate to study and model tsunamis as observed in nature. Therefore, it is of interest to generate tsunamis more appropriately in theoretical models and laboratory experiments. One way to achieve this is through the motion of a bottom disturbance on the seafloor.

In the literature of waves produced by a moving disturbance on the seafloor, various types of motions of the bottom disturbance are considered, including the bottom disturbance moving vertically, see, e.g., [13,14,15,16], or horizontally, see, e.g., [17,18,19], or along a slope boundary, see e.g., [20,21,22,23]. The effects of a moving surface disturbance (such as a ship) on the fluid field have been studied by, for example, [24,25,26,27], which also include the waves generated by the motion of the surface disturbance. In this study, we are focused on waves generated by a horizontally moving disturbance on a flat seafloor.

There have been a very limited number of laboratory experiments conducted to investigate wave generation by a horizontally moving disturbance on a flat seafloor. Ref. [18] conducted laboratory experiments, focusing on surface elevation generated by a horizontally moving semi-elliptical disturbance in a two-dimensional wave flume. It was found that the geometrical characteristics of the bottom disturbance significantly influence the generated wave amplitudes, while the acceleration of the bottom disturbance has less effect. Ref. [28] recorded the time series of surface elevation by the horizontally moving rectangular blocks within a three-dimensional wave tank, and three-dimensional effects were observed during the wave generation process.

Different theoretical models have been developed to study wave generation by horizontally moving disturbances on a flat seafloor, where the generated wave amplitude or energy are the parameters of interest. One-dimensional shallow-water-wave approximations were applied by, e.g., [17,29] to obtain the analytical solutions of wave generation by horizontally moving disturbances, where it was found that generated wave energy is proportional to the disturbance speed. Linear models were developed by [30,31], in which the fluid is governed by the Laplace equation subject to the linear boundary conditions. Within these studies, it is found that the geometric characteristics of the disturbance and the moving speed of the disturbance have a remarkable influence on the wave energy. Ref. [19] derived the analytical solutions for waves generated by a horizontally moving disturbance based on the shallow-water-wave equations. They showed that the enclosed area of the disturbance has more significant influence on the generated waves than the disturbance shape itself.

The effect of dispersion and nonlinearity on the generated waves by horizontally moving disturbances have also been theoretically investigated. Wave nonlinearity, defined as the ratio of wave height to wave flume depth, and wave dispersion, defined as wavelength to wave flume depth, change with varying disturbances. Based on the linear models, refs. [32,33] found that the dispersive effect was important for the case of intermediate and deep water depths. Ref. [34], by use of the linear and fully dispersive wave model, proposed that a peak energy exists during the wave generation process. Ref. [35] derived the GN equations to study the waves generated by a moving disturbance on the seafloor (disturbance define by an 8th-order polynomial). The results of the GN model were compared to the Boussinesq-class equations of [36,37], as well as to linear solutions. The study showed that the wave resistance, with a shifted peak value, fluctuated with the moving speed. The linear solutions predicted higher wave amplitudes compared to the other models, depending on the moving speed, and the Boussinesq models showed better agreement with the GN model in most cases. Ref. [38] developed a theoretical model based on the GN equations to repeat the same case of [31], and found that the proposed GN model’s results showed remarkable differences when compared to the linear solutions of [31] in certain cases, i.e., nonlinearity can play an important role. Ref. [39] using the computational fluid dynamics method observed that the laminar and turbulence models agreed well with the measured wave phases, but both laminar and turbulence results underestimated measured wave amplitudes. Ref. [40] used the moving particle explicit method and studied the effect of the nonlinearity and dispersion on the generated waves in an intermediate water depth.

It is also possible to generate periodic waves by the oscillatory motions of a disturbance on the seafloor. In this case, the moving disturbance acts as a wavemaker on the seafloor, forming oscillatory waves. The generated waves change form, depending on the shape of the bottom disturbance and the type of motion considered—see, e.g., [40,41,42].

In this study, we are attempting to investigate wave generation by a horizontally moving disturbance. The bottom disturbance generates long waves by a single motion or periodic waves by oscillatory motions. Nonlinearity and dispersion effects play important roles in the stability of the waves generated by the horizontally moving disturbances. The wave properties are studied by developing nonlinear, dispersive models to show how different parameters play a role in generating waves by the motion of a bottom disturbance, and how these wave properties, e.g., wave celerity and wave nonlinearity, are affected. Also, in this study, we consider two models, namely, an inviscid model based on the Level I GN equations and a Navier–Stokes (NS) model by solving the incompressible NS equations, and compare their results to assess how viscosity plays role and to what extend the proposed GN and NS models agree with each other. Studies on different shapes of the bottom disturbances and their interaction with generated waves, e.g., [43,44,45], are beyond the scope of this study, although this can be achieved through the proposed models; see, e.g., [46,47].

The GN equations and NS models are first introduced, followed by the description of the moving disturbance considered here. Then, the numerical setup and numerical solutions of the two models are introduced. Discussion on the nature of wave generation by a moving disturbance is given next by developing a model based on the linearized GN equations. In the results section, we first present the results for a single motion along with comparisons of the results of the two models with the available data. This is followed by a qualitative analysis of the pressure and velocity fields. A variety of parameters, including the geometric and kinematic characteristics of the disturbance, are considered, where the results of the time series of the surface elevation and variation of wave properties with different parameters are also investigated. The concluding remarks section closes the paper.

2. Theories

In this section, the two models developed to study wave generation by a horizontally moving disturbance on the seafloor are discussed. In these models, the fluid is assumed to be incompressible. A two-dimensional Cartesian right-handed coordinate system ($x_1$, $x_2$) is used, where $x_1$ points to the right and $x_2$ points vertically upwards, against the gravity acceleration. The coordinate system origin is located at the still-water level (SWL), shown in Figure 1. In all cases, the seafloor away from the moving disturbance is flat, though this is not required in general.

2.1. The GN Equations

Refs. [48,49] proposed the GN theory for wave propagation in an incompressible, inviscid fluid sheet. The GN equations describe nonlinear fluid motion, while they satisfy the conservation of mass and nonlinear boundary conditions exactly and postulate the conservation of momentum in an integrated manner. No restriction is needed about the flow to be irrotational. The only assumption about the fluid kinematics is regarding the vertical velocity distribution. In the Level I GN equations, the vertical velocity is distributed linearly in the vertical direction, resulting in invariant horizontal velocity over the vertical column in this inviscid fluid sheet.

Ref. [35] provided the Level I GN equations in a compact form, given as

$${\eta }_{,t}+{\left\{(h+\eta -\alpha )u_1\right\}}_{,x_1}={\alpha }_{,t},$$

(1)

$${\dot{u}}_1+g{\eta }_{,x_1}+{\displaystyle \frac{{\widehat{p}}_{,x_1}}{\rho }}=-{\displaystyle \frac{1}{6}}\{{[2\eta +\alpha ]}_{,x_1}\ddot{\alpha }+{[4\eta -\alpha ]}_{,x_1}\ddot{\eta }+(h+\eta -\alpha ){[\ddot{\alpha }+2\ddot{\eta }]}_{,x_1}\},$$

(2)

$$u_2(x_1,x_2,t)=\dot{\alpha }+{\displaystyle \frac{x_2+h-\alpha }{h+\eta -\alpha }}(\dot{\eta }-\dot{\alpha }),$$

(3)

$$\overline{p}(x_1,t)={\displaystyle \frac{\rho }{2}}(h+\eta -\alpha )(\ddot{\alpha }+\ddot{\eta }+2g)+\widehat{p},$$

(4)

where $\eta (x_1,t)$ is the surface elevation on the top boundary, measured from the SWL, $\alpha (x_1,t)$ is the deformation on the bottom boundary, h is the constant water depth, $u_1(x_1,x_2,t)$ and $u_2(x_1,x_2,t)$ are the horizontal and vertical fluid particle velocities, respectively, g is the gravitational acceleration, $\rho$ is the fluid density, $\widehat{p}(x_1,t)$ is the fluid pressure on the top surface, taken as $\widehat{p}=0$ without losing generality and $\overline{p}(x_1,t)$ is the pressure on the bottom boundary. Subscripts after commas denote partial differential derivatives with respect to corresponding quantities. For arbitrary variable $\beta$, $\dot{\beta }$ represents the first-order total (material) derivatives and $\ddot{\beta }$ represents the second-order total derivatives.

Following [50], the fluid pressure distribution in the domain can be obtained as

$$p(x_1,x_2,t)={\displaystyle \frac{\rho }{2}}(h+\eta -\alpha )(\ddot{\alpha }+\ddot{\eta }+2g)-\rho g(x_2+h-\alpha )-{\displaystyle \frac{\rho }{2}}{(x_2+h-\alpha )}^2{\displaystyle \frac{\ddot{\alpha }+\ddot{\eta }}{h+\eta -\alpha }},$$

(5)

where $p=0$ at $x_2=\eta$ (assuming that the top pressure is atmosphere) and $p=\overline{p}$ at $x_2=-h+\alpha$.

Expanding the material derivatives, Equation (2) is expressed as

$$\begin{array}{cc}\hfill & {u_1}_{,t}+u_1{u_1}_{,x_1}+g{\eta }_{,x_1}=\hfill \\ \hfill & -{\eta }_{,x_1}[{\alpha }_{,tt}+{u_1}_{,t}{\alpha }_{,x_1}+2u_1{\alpha }_{,x_{1}t}+u_1{u_1}_{,x_1}{\alpha }_{,x_1}+u_1^2{\alpha }_{,x_1{x}_1}]\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}(h+\eta -\alpha )[{\alpha }_{,x_{1}tt}+{u_1}_{,x_{1}t}{\alpha }_{,x_1}+{u_1}_{,t}{\alpha }_{,x_1{x}_1}+2{u_1}_{,x_1}{\alpha }_{,x_{1}t}\hfill \\ \hfill & +2u_1{\alpha }_{,x_1{x}_{1}t}+{u_1^2}_{,x_1}{\alpha }_{,x_1}+u_1{u_1}_{,x_1{x}_1}{\alpha }_{,x_1}+3u_1{u}_{,x_1}{\alpha }_{,x_1{x}_1}+u_1^2{\alpha }_{,x_1{x}_1{x}_1}]\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}(h+\eta -\alpha )({u_1^2}_{,x_1}-{u_1}_{,x_{1}t}-u_1{u_1}_{,x_1{x}_1})\left[{(2\eta -\alpha )}_{,x_1}\right]\hfill \\ \hfill & -{\displaystyle \frac{1}{3}}{(h+\eta -\alpha )}^2({u_1}_{,x_1}{u_1}_{,x_1{x}_1}-{u_1}_{,x_1{x}_{1}t}-u_1{u_1}_{,x_1{x}_1{x}_1}).\hfill \end{array}$$

(6)

The term $\alpha$ and its spatial and time partial derivatives in Equations (1) and (6) are defined by the motion of the bottom disturbance, which are explicitly given and are discussed in the following sections. In Equations (1) and (2), the unknowns are $\eta$ and $u_1$; thus, the number of unknowns (2) is equal to the number of equations to close the system. Various previous studies, e.g., [51,52,53,54], have shown the accuracy and efficiency of the Level I GN equations, Equations (1) to (4), to analyze nonlinear wave mechanics and wave–body interactions in shallow waters. For the higher-level GN model, its vertical velocity distribution is in the form of more complex formulas, such as high-order polynomial functions; see, e.g., [55,56,57,58].

2.2. The Navier–Stokes (NS) Equations

In the NS model, the fluid is assumed to be incompressible and homogeneous, and the laminar fluid motion is governed by the incompressible NS equations, described by

$$u_{j,j}=0,j=1,2,$$

(7)

$$u_{i,t}+u_j{u}_{i,j}=-{\displaystyle \frac{1}{\rho }}{p}_{,i}+\nu u_{i,jj}+g_i,i,j=1,2,$$

(8)

where p represents the fluid pressure, $g_i=(0,-g)$ and $\nu$ represents the fluid kinematic viscosity. Einstein’s summation convention is used in the indicial notation for the repeated indices. In this study, we assume that the turbulence effect is negligible and, thus, the laminar model is used, which was successfully applied for wave interaction with submerged structures by [59,60].

The unknowns in Equations (7) and (8), i.e., p, $u_1$ and $u_2$, are simultaneously solved for the air and water phases, and the free surface between the two phases is captured by applying the volume of the fluid method (VoF); see, e.g., [61,62,63].

3. The Moving Disturbance

In general, the bottom disturbance can be of an arbitrary shape and move in a single or multiple oscillatory motions. In this study, we are focused on a continuously and horizontally moving disturbance on the seafloor. The disturbance on the seabed, $\alpha$, is given by an arbitrary continuous formula of $x_1$ and t:

$$\alpha =f(x_1,t),$$

(9)

where $f(x_1,t)$ can be any type of continuous function in general. In this study, we use the following function to describe the bottom disturbance:

$$f(x_1,t)=A_0\times R(x_1,t)\times D(x_1,t),$$

(10)

where $A_0$ is the geometric amplitude for the bottom disturbance as shown in Figure 1, and

$$D(x_1,t)=sec{h}^2(x_1-x_0),$$

(11)

where $x_0=x_0\left(t\right)$ is the instantaneous position of the center of the moving disturbance. $R(x_1,t)$ is the ramp function to create smooth edges, described as

$$R(x_1,t)=e^{{\displaystyle \frac{-{(x_1-x_0)}^2}{{\sigma }^2}}},$$

(12)

where $\sigma ={\displaystyle \frac{L_s}{4}}$ is constant, defining the effective length of $R(x_1,t)$, and $L_s$ is the disturbance length, shown in Figure 1.

We consider both a single motion and oscillatory motion of the bottom disturbance. For the single motion, to eliminate the effect of the acceleration of the bottom disturbance, we assume that the bottom disturbance starts moving constantly at $t=0$, and the disturbance’s speed progresses to zero at the end of motion. Hence, $x_0\left(t\right)$ is described by

$$x_0\left(t\right)=\left\{\begin{array}{c}{X}_0+V_{s}t,t\le T_0,\\ X_0+L_h,t>T_0,\end{array}\right.$$

(13)

where $X_0$ is the initial coordinate of the disturbance center, $L_h$ is the moving distance of the bottom disturbance, $T_0$ is the movement duration and $V_s$ is the speed. The motion distance $L_h$ is determined by

$$L_h=V_s{T}_0.$$

(14)

Physically, the velocity of the bottom disturbance would vary gradually at the beginning and end of the motion, and that would alter the form of the waves generated—see e.g., [18,64].

For the oscillatory motion of the bottom disturbance, $x_0\left(t\right)$ is described by

$$x_0\left(t\right)=X_0+A_h\times s\left(t\right)\times sin\left(\omega t\right),$$

(15)

where $A_h$ is the oscillation amplitude, $\omega$ is the frequency of the oscillation and $s\left(t\right)$ is the smoothing function used to control the oscillations, which is described by

$$s\left(t\right)=1-e^{-qt},$$

(16)

where q is a constant used to define a gradually increasing oscillation (ramp). It is found that $q=0.1$ provides a stable solution while keeping the computational cost reasonable, which is similar to that used by [65].

4. Numerical Setup and Solution

In this section, the numerical setup and solution of the GN and NS models are discussed. All variables and results in this manuscript are given in a non-dimensional form by g, h and $\rho$ from now on, such that

$$\begin{array}{cc}\hfill & x_1^{\prime }={\displaystyle \frac{x_1}{h}},{\eta }^{\prime }={\displaystyle \frac{\eta }{h}},{\alpha }^{\prime }={\displaystyle \frac{\alpha }{h}},\hfill \\ \hfill & u_1^{\prime }={\displaystyle \frac{u_1}{\sqrt{gh}}},t^{\prime }=t\sqrt{{\displaystyle \frac{g}{h}}},p^{\prime }={\displaystyle \frac{p}{\rho gh}}.\hfill \end{array}$$

(17)

The superscripts (′) are removed from all dimensionless variables for simplicity.

4.1. The GN Modeling

Substituting the form of the bottom disturbance considered in this study, given by (10), the spatial and time partial derivatives of function $\alpha$ in Equations (1) and (6) read as

$$\begin{array}{c}\hfill {\alpha }_{,t}=A_0\times (R_{,t}D+R{D}_{,t}),\end{array}$$

(18a)

$$\begin{array}{c}\hfill {\alpha }_{,x_1}=A_0\times (R_{,x_1}D+R{D}_{,x_1}),\end{array}$$

(18b)

$$\begin{array}{c}\hfill {\alpha }_{,tt}=A_0\times (R_{,tt}D+2R_{,t}{D}_{,t}+R{D}_{,tt})\end{array}$$

(18c)

$$\begin{array}{c}\hfill {\alpha }_{,x_{1}t}=A_0\times (R_{,x_{1}t}D+R_{,t}{D}_{,x_1}+R_{,x_1}{D}_{,t}+R{D}_{,x_{1}t}),\end{array}$$

(18d)

$$\begin{array}{c}\hfill {\alpha }_{,x_1{x}_{1}t}=A_0\times (R_{,x_1{x}_{1}t}D+2R_{,x_{1}t}{D}_{,x_1}+2R_{,x_1}{D}_{,x_{1}t}+R_{,t}{D}_{,x_1{x}_1}+R_{,x_1{x}_1}{D}_{,t}+R{D}_{,x_1{x}_{1}t}),\end{array}$$

(18e)

$$\begin{array}{c}\hfill {\alpha }_{,x_{1}tt}=A_0\times (R_{,x_{1}tt}D+2R_{,x_{1}t}{D}_{,t}+2R_{,t}{D}_{,x_{1}t}+R_{,tt}{D}_{,x_1}+R_{,x_1}{D}_{,tt}+R{D}_{,x_{1}tt}),\end{array}$$

(18f)

$$\begin{array}{c}\hfill {\alpha }_{,x_1{x}_1}=A_0\times (R_{,x_1{x}_1}D+2R_{,x_1}{D}_{,x_1}+R{D}_{,x_1{x}_1}),\end{array}$$

(18g)

$$\begin{array}{c}\hfill {\alpha }_{,x_1{x}_1{x}_1}=A_0\times (R_{,x_1{x}_1{x}_1}D+3R_{,x_1{x}_1}{D}_{,x_1}+3R_{,x_1}{D}_{,x_1{x}_1}+R{D}_{,x_1{x}_1{x}_1}).\end{array}$$

(18h)

Relations (18a)–(18h) are substituted into Equations (1) and (6).

Wave absorbers are applied on the two side boundaries of the wave tank by applying the Orlanski condition (see, e.g., [66]), given by

$${\mathsf{\Omega }}_t+c{\mathsf{\Omega }}_{x_1}=0,$$

(19)

where c is the linear wave celerity and $\mathsf{\Omega }$ represents $\eta$ or $u_1$. The tank length is also long enough to minimize wave reflections into the domain. Although this is an approximate condition, it has proven to be quite accurate in previous studies as long as the tank length is long; see, e.g., [24].

A second-order spatial central-difference method is used to obtain the numerical solution of the governing equations. The unknowns of the GN model in Equations (1) and (6) are $\eta (x_1,t)$ and $u_1(x_1,t)$. Within the finite-difference method, all continuous quantities, $\eta (x_1,t)$ and $u_1(x_1,t)$, are represented by the discretized values $\eta (i,n)$ and $u_1(i,n)$, respectively, where i and n are the $i^{th}$ and $n^{th}$ spatial mesh points and time steps, respectively. For time marching, the Modified Euler Method (MEM) is applied.

The two-step MEM calculation is first performed to obtain $\eta$ and $u_1$. Then, $x_0$, R and D are updated to determine $\alpha$. The process is repeated in each time step, and it continues for $n_{max}$ times. The workflow of the numerical procedure of the GN model is shown in Figure 2.

4.2. The NS Modeling

A two-dimensional, rectangular computational tank, similar to that shown in Figure 1, is created to study the generation of surface waves by bottom disturbances. The water depth in the NS computational tank is invariant, and $h=1$ m for all cases considered here. Within the tank, the motion of the fluid is governed by the NS equations. Appropriate boundary conditions are enforced, including the no-slip wall condition on the seafloor and atmospheric pressure on the top of the domain. For both the left and right boundaries of the tank, a no-slip wall boundary condition is applied.

The NS equations are discretized and solved within the framework of the finite-volume method. The computational region is divided into a series of non-overlapping control volumes. The the discrete values located at the center of the control volume are used to approximate continuous physical quantities.

To resolve the incompressible NS equations, the PIMPLE algorithm is used for the pressure–velocity coupling, which is a merged algorithm of Pressure Implicit with Splitting of Operators and Semi-Implicit Method for Pressure-Linked Equations; see, e.g., [67].

The dynamic mesh technique is applied to update the mesh at each time step for the motion of the bottom disturbance. Mesh deformations of the bottom boundary are prescribed by Equations (13) and (15) for the single and oscillatory motions, respectively. A mesh morphing technique is used to determine the new mesh of the entire domain. Mesh point distortion on the bottom boundary is interpolated into the mesh cell vertices inside the domain; see, e.g., [68,69] for more details.

The flowchart of the numerical procedure of the NS model is presented in Figure 3. The velocity and pressure fields are obtained in each time step. The surface elevation between the air and water phases are obtained by the use of the VoF method. The oscillations of the bottom disturbance are updated and the mesh deformations on the bottom are found by the mesh morphing approach. Thus, new mesh points are updated over the entire domain.

A mesh convergence analysis is performed for the case of a single motion of the bottom disturbance where $L_s=4$, $A_0=0.05$, $V_s=0.5$, $L_h=10$, $T_0=20$ and the tank length is 70. Four numerical wave gauges are placed in the tank, whose locations are shown in Figure 1. For the case of the single motion, the locations of the four gauges shown in Figure 1 are fixed at $L_1=1.2L_h$, $L_2=0.8L_h$ and $L_3=L_h$.

Table 1. Mesh configurations used for the mesh convergence study. The length and height of the computational domain are 70 and 0.5, respectively.

Configurations	$\mathbf{\Delta }{\mathit{x}}_1$	$\mathbf{\Delta }{\mathit{x}}_2$ (Above and Below the Free Surface)	CPU Time (Hours)
Mesh 1	$L_s$/50	$A_0$/10	0.91
Mesh 2	$L_s$/100	$A_0$/10	3.22
Mesh 3	$L_s$/200	$A_0$/10	12.74
Mesh 4	$L_s$/100	$A_0$/5	3.69
Mesh 5	$L_s$/100	$A_0$/20	3.68

presents different mesh configurations applied in the mesh convergence study. The grid size in the horizontal direction, $\mathsf{\Delta }{x}_1$, is uniform while mesh size around the free surface in the vertical direction, $\mathsf{\Delta }{x}_2$, is particularly refined. All computations are performed for the same simulation time ($t=75$) in OpenFOAM, and by a workstation of twenty cores, 2.20 GHz Intel Xeon Silver 4114 and 32 GB memory.

Time step ($\mathsf{\Delta }t$) is adjusted according to the maximum Courant number, $C{o}_{max}$. The Courant number is given by $Co=|u_i|\mathsf{\Delta }t/\mathsf{\Delta }{x}_i$, where the subscripts $i=1$ and 2 refer to the components in the horizontal and vertical directions, respectively. $|u_i|$ is the amplitude of the velocity vector, and the spatial grid size is $\mathsf{\Delta }{x}_i$.

The time series of surface elevation at Gauges GII and GIV are shown in Figure 4. It is found that the mesh configuration Mesh 2 and $C{o}_{max}=0.1$ is converged particularly for the first wave peaks and it is appropriate to perform the computations. Thus, this configuration is applied for the case of the bottom disturbance with single motion in the following sections. We perform the same convergence study for the oscillatory motion and the converged configuration is $\mathsf{\Delta }{x}_1=L_s/160$, $\mathsf{\Delta }{x}_2=A_0/10$ and $C{o}_{max}=0.1$.

Similar mesh convergence analysis is performed in the GN model, and mesh configuration $\mathsf{\Delta }{x}_1=0.05$ and $\mathsf{\Delta }t=0.025$ is found to be suitable among all configurations considered in the study. The computations of the GN equations for the cases considered here are concluded within minutes on a laptop with a single CPU of 1.60 $\mathrm{GHz}$ Intel(R) Core(TM) i5-8265U CPU and 16 GB memory.

5. Wave Generation by the Motion of a Bottom Disturbance

In this section, we discuss the mechanism for wave generation by a moving disturbance on the seafloor. For this purpose, we will first obtain the linearized form of the GN equations to access the generation of linear waves.

Following the approach given in [35,65] by neglecting nonlinear terms in Equations (1) and (2), after expanding all material derivative terms, the linear GN (LGN) equations are obtained as

$${\eta }_{,t}+{u_1}_{,x_1}={\alpha }_{,t},$$

(20)

$${u_1}_{,t}+{\eta }_{,x_1}=-{\displaystyle \frac{1}{6}}{\alpha }_{,tt{x}_1}-{\displaystyle \frac{1}{3}}{\eta }_{,tt{x}_1}.$$

(21)

We apply the LGN equations to study how the motion of a disturbance on the seafloor results in the generation of surface waves.

Rearranging Equation (20) gives

$$\begin{array}{c}\hfill {\eta }_{,t}={\alpha }_{,t}-{u_1}_{,x_1}.\end{array}$$

(22)

$\eta$ can be numerically obtained by use of the finite-difference method, determined by

$$\begin{array}{c}\hfill {\eta }^{n+1}={\eta }^n+({\alpha }_{,t}-{u_1}_{,x_1})\mathsf{\Delta }t,\end{array}$$

(23)

where ${\eta }^n$ and ${\eta }^{n+1}$ are the surface elevation in the current and next time steps. Equation (23) shows explicitly that ${\eta }^{n+1}$ is a function of $\alpha$, i.e., that the deformation of the surface elevation is caused by the moving disturbance on the seafloor.

Equation (21) is the linear form of Equation (2). Taking partial differentiation on both sides of Equation (20) and rearranging the relation gives

$${\eta }_{,tt{x}_1}={\alpha }_{,tt{x}_1}-u_{1,x_1{x}_{1}t}.$$

(24)

Substitute Equation (24) into Equation (21), we write

$${u_1}_{,t}-{\displaystyle \frac{1}{3}}{u}_{1,x_1{x}_{1}t}=-{\eta }_{,x_1}-{\displaystyle \frac{1}{2}}{\alpha }_{,tt{x}_1}.$$

(25)

$u_1$ can be obtained by numerically solving (25) by use of the finite-difference method.

Expanding all material derivative terms in Equation (3) and neglecting the nonlinear terms, the linear vertical velocity reads as

$$\begin{array}{c}\hfill u_2(x_1,x_2,t)={\alpha }_{,t}+{\displaystyle \frac{x_2+1-\alpha }{1+\eta -\alpha }}({\eta }_{,t}-{\alpha }_{,t}).\end{array}$$

(26)

The vertical velocity on the seafloor and the free surface can be obtained by substituting $x_2=-1+\alpha$ and $x_2=\eta$ into Equation (26), respectively. $u_2$ is then described by

$$u_2(x_1,-1+\alpha ,t)={\alpha }_{,t},$$

(27)

$$u_2(x_1,\eta ,t)={\eta }_{,t},$$

(28)

where ${\alpha }_{,t}$ and ${\eta }_{,t}$ are the exact kinematic boundary conditions on the bottom and the free surface, respectively, for the LGN model. From Equation (28), $\eta =0$ if $u_2(x_1,\eta ,t)=0$. Equation (28) confirms explicitly that the wave surface elevation is deformed due to the vertical momentum around the free surface.

We recall that, in the Level I GN equations in this study, $u_2$ is distributed linearly in the vertical direction. Thus, from Equations (27) and (28), i.e., velocities at the bottom and top surfaces of the fluid sheet, we can define the distribution of $u_2$ along the water column. Figure 5 shows the vertical distributions of $u_2$ at the trailing edge of the bottom disturbance, center and the leading edge of the bottom disturbance by the GN and LGN models for $L_s=4$, $A_0=0.1$, $V_s=0.5$ and $L_h=10$. Shown in Figure 5, we observe different slopes at the trailing edge, center and leading edge of the bottom disturbance in the vertical direction at different times.

Shown in Figure 5, we find that $u_2$ at the leading edge of the moving disturbance is positive while $u_2$ at its trailing edge is negative. It is also observed that $u_2$ at the center of the moving disturbance is near zero at all times. That is, the leading edge of the disturbance pushes the fluid upwards, while the trailing edge causes a downward motion.

Shown in Figure 5, $u_2$ on the seafloor obtained by the LGN model is quite close to the nonlinear GN model. We also find that the GN and LGN models predict similar results at the leading edge of the moving disturbance in this case. The largest differences are observed where nonlinearity becomes more dominant at the trailing edge and center of the moving disturbance, particularly at $t\ge 10$, where nonlinearity becomes more significant.

The snapshots of $\alpha$, $u_2$ on the seafloor, $u_2$ on the free surface and $\eta$ obtained by the GN and LGN models at $t=1$, $t=5$, $t=10$ and $t=20$ are shown in Figure 6 for $L_s=4$, $A_0=0.1$, $V_s=0.5$ and $L_h=10$. We find that $u_2$ on the seafloor obtained by the LGN model is very close to that of the GN model. Shown in Figure 6, a leading peak is formed on the free surface, followed by a series of tail waves. GN and LGN models predict very close agreement on the leading-wave peaks because the nonlinearity for $u_2$ is less significant above the leading edge of the disturbance, shown in Figure 5 as well.

6. Results: Single Motion of the Bottom Disturbance

Results obtained by the GN and NS models are presented in this section, starting with comparisons with available data and the results of pressure and velocity fields for the single-motion condition.

The effect of the disturbance amplitude ($A_0$) and disturbance speed ($V_s$) on the generated waves are presented and discussed in Section 6.3 and Section 6.4. Our analysis indicates that disturbance length ($L_s$) and the disturbance motion range ($L_h$) have little effect on the results and, therefore, they are not presented in this study.

6.1. Comparisons with the Available Data

Ref. [17] used shallow-water approximations to study wave generation by a bottom disturbance of single motion. This case is used here for comparison purpose. In this case, the geometric shape of the bottom disturbance is described by

$$\alpha (x_1,t)=\left\{\begin{array}{c}{\displaystyle \frac{A_0}{2}}[1-cos\left({\displaystyle \frac{2\pi }{L_s}}(x_1-x_s)\right)],x_s\le x_1\le x_e,\\ 0,x_1<x_s,x_1>x_e.\end{array}\right.$$

(29)

This case is studied here by use of the GN and NS models and for the following conditions: $L_s=10$, $A_0=0.01$, $V_s=0.639$ and $T_0=25.06$. The snapshots of surface elevation at different times obtained by the two models and the results of [17] are shown in Figure 7. Figure 7 shows that two waves are generated due to the motion of the bottom disturbance: one wave with a distinguished peak and trough propagates in the same direction as the moving bottom, although a single depression propagates in the opposite direction. We find that the following tail waves are formed by the nonlinear models, which are omitted in the linearized model in [17]. Those tail waves due to nonlinearity modulate the leading waves and affect wave evolution: the peak and trough obtained by the two models propagate slightly more slowly than that in [17]. Overall, quite good agreement can be found between the two models and the results in [17]. The GN and NS models show closer agreement with each other.

6.2. Snapshots of the Velocity and Pressure Fields

In this section, the evolution of the pressure and velocity fields are presented and discussed for the case of $L_s=4$, $A_0=0.1$, $V_s=0.5$ and $L_h=10$. The snapshots of the pressure and velocity vector fields at $t=1$, $t=2$, $t=5$ and $t=10$ obtained by the NS and GN models are shown in Figure 8 and Figure 9, respectively.

Shown in the figures, the leading edge of the bottom disturbance results in an upward motion of the fluid, while its trailing edge causes a downward motion. Shown in Figure 8 and Figure 9, wave surface elevation is found to be slightly disturbed at $t\le 2$, while more notable disturbances are formed at $t\ge 5$. This is because the velocity vectors around the free surface are smaller at $t\le 2$ and larger at $t\ge 5$. Overall, we find that the GN model shows comparable agreement with the NS model.

The evolution of the generated wave surface elevation by a moving disturbance at $t=1$, $t=2$, $t=5$, $t=10$, $t=15$ and $t=20$ is shown in Figure 10 for the case of $L_s=4$, $A_0=0.1$, $V_s=0.5$ and $L_h=10$. It is observed that a leading-wave peak followed by a series of tail waves are generated during the process. Shown in the figure, the leading-wave peak propagates faster than the moving disturbance at $t\ge 10$. The GN models agrees well with the NS model for the leading wave at $t\ge 5$.

6.3. Effect of Disturbance Amplitude

In this subsection, we study how the wave field is affected by the change in disturbance amplitude. The time series of surface elevation for $A_0=0.05$, $A_0=0.10$ and $A_0=0.15$ and $L_s=4$, $V_s=0.5$ and $L_h=10$ are shown in Figure 11. The surface elevation is recorded at Gauges GI, GII, GIII and GIV, and their positions are shown in Figure 1. $L_1=1.2L_h$, $L_2=0.8L_h$ and $L_3=L_h$. The maximum disturbance amplitude $A_0=0.15$ is used for long waves in the study, but this is not the limitation of the present models in general.

Figure 11 shows that a main leading wave with a series of oscillatory tail waves are formed. For different $A_0$, the wave phases generated are similar, while the generated wave amplitudes are different. The largest differences between the two models are observed in the tail dispersive waves. Overall, the two models predict close results for leading-wave peaks in all cases.

Properties of the leading-wave are particularly analyzed in the study, and we define A as the wave amplitude of the leading wave. A is obtained by determining the first peak value from the time series of the wave surface elevation in different wave gauges. The variation of A with different $A_0$ is shown in Figure 12. We find that A increases nonlinearly with $A_0$ in all cases. We observe that A generally decreases with longer propagating distance, but this is less remarkable for $A_0=0.05$.

6.4. Effect of the Motion Speed

The effect of the disturbance speed on the wave field is discussed in this section. The time series of surface elevation for $V_s=0.4$, $V_s=0.6$ and $V_s=0.8$ are shown in Figure 13. A leading wave is produced, followed by a series of tail waves. We find that the wave free surface is very different for different $V_s$. It is observed that the tail waves upwave decrease over time for $V_s=0.6$ and $V_s=0.8$. The tail waves downwave are dispersive and behave like periodic waves for $V_s=0.8$. Overall, the two models show close results, particularly for the leading waves.

To determine whether the leading waves generated are propagating at critical or supercritical speeds, like solitary waves, at this point, we assess the propagation speed of the leading-wave by tracking leading wave peaks between wave gauges—see, e.g., [51] on the speed of solitons. The analytical relation for the propagation speed of a solitary wave, U, is given by (see, e.g., [70])

$$U=\sqrt{A+1},$$

(30)

in the Level I GN theory. Thus, $U\ge 1$ is satisfied for a soliton; see, e.g., [71].

In this study, the speed of the leading-wave is determined from the first peaks of the time series of wave surface elevation at different gauges. We use $U_1$, $U_2$ and $U_3$ to represent the leading-wave speed at areas between GI and GII, GII and GIII and GIII and GIV, respectively.

The variation of $U_1$, $U_2$ and $U_3$ with $V_s$ is given in

Table 2. Variation of the propagation speed of the leading wave with $V_s$.

${\mathit{V}}_{\mathit{s}}$	${\mathit{U}}_1$ (Between GI and GII)	${\mathit{U}}_2$ (Between GII and GIII)	${\mathit{U}}_3$ (Between GIII and GIV)
0.4	1.025	1.023	1.005
0.5	1.031	1.015	1.023
0.6	1.025	1.005	1.005
0.7	0.972	0.995	0.998
0.8	0.969	0.983	0.995

. We find that the leading waves propagate with supercritical speed, like solitons, in all cases when $V_s=0.4$, $V_s=0.5$ and $V_s=0.6$. For the case of $V_s=0.7$ and $V_s=0.8$, we find that the wave propagation speeds are below critical (less than one). This illustrates that, when the bottom disturbance moves with $V_s\ge 0.7$, the generated leading waves propagate differently from solitary waves.

7. Results: Oscillatory Motion of the Bottom Disturbance

In this section, waves generated by an oscillatory motion of the bottom disturbance, described by Equation (15), are investigated by using the GN and NS models. Results, including snapshots of the velocity and pressure fields and wave surface elevation, are followed by a parametric study considering a range of disturbance amplitudes and oscillation amplitudes.

Our assessment shows that the disturbance length ($L_s$) and oscillation period (T) play negligible roles on the generated waves and, thus, are not included in the analysis. Wave gauge locations are fixed at $L_1=L_2=L_3=4$.

7.1. Snapshots of the Velocity and Pressure Fields

In this section, we show the evolution of the pressure and velocity fields generated by the oscillatory motion of the bottom disturbance over one period for the case of $L_s=4$, $A_0=0.1$, $A_h=0.8$ and $T=6$. The snapshots of the pressure and velocity fields from $t=72$ to $t=77.25$, obtained by the NS and GN models, are shown in Figure 14 and Figure 15, respectively. The corresponding time series of the oscillatory center of the bottom disturbance is specified as $x_0=X_0$ at $t=72$ and $t=75$, approximately.

Shown in Figure 14 and Figure 15, the bottom disturbance causes the surrounding fluid to undergo oscillatory motions in the vertical direction. It is observed that the oscillatory bottom disturbance modifies the velocity vectors over the free surface significantly, and this causes deformations around the free surface.

A snapshot of the wave profile obtained by the two models at $t=125.3$ for $L_s=4$, $A_0=0.1$, $A_h=0.8$ and $T=6$ is shown in Figure 16. In this figure, we are mainly interested in waves approximately within the region $10\le x_1-X_0\le 20$, within which the waves are fully developed and outside of the ramp and the wavefront zones. It is observed that the peaks of the GN waves are very close to the Airy wave theory, while the NS waves are slightly smaller, possibly due to the viscous effects. The linear waves propagate slightly more slowly than the GN waves, but their differences are not significant. The NS waves also advance a little more quickly than the linear waves and the differences are more noticeable. We note that part of the differences observed between the NS waves and others might be due to the effect of the developing wavefront and the transient waves; see, e.g., [72,73] for more details. The use of a substantially longer domain and significantly longer duration for the wave generation would be ideal; however, this is not practical with the NS model due to the need for substantially larger computational power and time required. For further discussion on establishing an efficient and practical numerical wave tank for NS equation, refer to, e.g., [74].

The vertical distribution of $u_1$, under the marked wave peak and trough, is shown in Figure 17, and distribution of the vertical velocity along the water column under the position of $\eta =0$, approximately between the marked wave peak and trough, is also presented in the same figure. We find that significant differences are observed for the $u_1$ distribution using the GN equations, Airy wave and NS equations. Compared to the Airy wave and NS approaches, $u_1$ calculated by the GN is constant in the vertical direction. In the NS model, $u_1$ is prescribed as zero on the bottom of the seafloor due to the no-slip boundary condition on the flat-bottom seafloor, i.e., it is stationary away from the disturbance. For the NS model and Airy wave solutions, $u_1$ decreases nonlinearly with the increasing water depth. In both cases, the NS model shows smaller $u_1$ because it predicts smaller wave peaks and troughs, as shown in Figure 16. $u_2$ is distributed linearly between the seabed and free surface for the proposed GN model, as expected. The other two models also show a nearly linear distribution in the vertical direction for $u_2$.

7.2. Effect of Disturbance Amplitude

In this subsection, we study how the waves generated by the oscillation of the bottom disturbance are affected by the disturbance amplitude. The time series of surface elevation for $A_0=0.1$, $A_0=0.2$ and $A_0=0.3$ and $L_s=4$$A_h=0.8$ and $T=6$ are shown in Figure 18. The generated waves propagate stably in the domain with no change in form or shape as a result of the balance between nonlinearity and dispersion. In all cases, the generated wave period exhibits less than 1 % deviation from the oscillation period of the bottom disturbance. We find that the deformations of wave surface elevation vary significantly with $A_0$. The differences between the two models are less remarkable for $A_0\le 0.2$, while more significant differences are found in $A_0=0.3$.

${\eta }_{max}$ and ${\eta }_{min}$ are determined from the wave peaks and troughs of the wave surface-elevation time series, respectively, in which ${\eta }_{max}$ and ${\eta }_{min}$ are averaged over three cycles. The variation of ${\eta }_{max}$, ${\eta }_{min}$ and wave height H ($H={\eta }_{max}-{\eta }_{min}$) with $A_0$ for Gauge GI are shown in Figure 19. ${\eta }_{max}$ and ${\eta }_{min}$ vary almost linearly with $A_0$ for the GN model, while they vary nonlinearly for the NS model. Shown in the figure, larger waves, i.e., larger wave peaks and troughs, are generated by the oscillatory motion of the bottom disturbance with a larger $A_0$. Thus, $A_0$ has a positive influence on H, i.e., wave nonlinearity increases with a larger $A_0$. Very close agreement between the two models is found for ${\eta }_{min}$ at $A_0\le 0.2$. For $A_0=0.3$, the NS model shows similar peaks but much smaller troughs than the GN model, i.e., waves produced by the NS model have relatively flat troughs, resulting in slightly smaller wave heights.

7.3. Effect of Oscillation Amplitude

The waves generated vary with the oscillation amplitude, and this is studied in this section. The time histories of wave elevation for $A_h=0.6$, $A_h=0.8$ and $A_h=1.0$ and $L_s=4$, $A_0=0.1$ and $T=6$ are shown in Figure 20. Stable, periodic waves are generated by the oscillatory motion of the bottom disturbance. Slight differences between the two models are found for different $A_h$. The two models agree very well on surface-elevation upwave calculations.

The variation of ${\eta }_{max}$, ${\eta }_{min}$ and H with $A_h$ for Gauge GI is shown in Figure 21. Shown in the figure, ${\eta }_{max}$ and ${\eta }_{min}$ vary nonlinearly with $A_h$ in all cases. It is observed that $A_h$ has a positive effect on H, i.e., H increases with increasing $A_h$ and this is nonlinear. We find that the ${\eta }_{max}$ and the absolute value of ${\eta }_{min}$ are increasing for larger $A_h$, i.e., waves generated have larger wave heights when the bottom disturbance oscillates with longer distance.

8. Concluding Remarks

Wave generation by a horizontally moving disturbance on an otherwise flat seafloor is studied by developing two theoretical models, namely, the GN and NS models. Both the single and oscillatory horizontal motions of the bottom disturbance are considered.

The linearized GN equations are proposed to better understand the mechanism of wave generation by a moving disturbance on the seafloor. We find that the nonlinear terms significantly affect the vertical velocity distributions above the trailing edge and center of the bottom disturbance.

A leading-wave followed by a series of tail waves are generated for a single motion, while periodic linear or nonlinear waves are generated by an oscillatory motion. Wave profiles generated by the oscillatory disturbance are found to be stable due to the balance between the wave nonlinearity and dispersion. Compared to the GN model, the NS model predicts larger leading-wave amplitude for a single motion but smaller wave amplitudes for an oscillatory motion.

A parametric study is performed for the single motion and oscillatory motion of the bottom disturbance, causing different geometries and motions of the bottom disturbance. The effect of the involved parameters are investigated. For the single-motion case, the leading-wave amplitudes vary nonlinearly with the involved parameters. The GN equations predict close results to the NS model in tail waves in most cases. The NS model predicts slightly larger leading wave peaks than the GN model in most cases considered in this study. The speed of the moving disturbance, $V_s$, significantly affects the generated leading waves. The generated leading-wave speed is smaller than that of solitons for $V_s\ge 0.7$.

For the oscillatory motion case, generated wave amplitudes vary linearly with the geometrical amplitude of the bottom disturbance, $A_0$, but vary nonlinearly with oscillation amplitude ($A_h$). The wave period is almost identical to the oscillation period of the bottom disturbance in all cases. We find that $A_0$ has an essential effect on the nonlinearity of the waves. Compared to the GN model, the NS model shows smaller downwave wave heights.

Overall, the agreement between the inviscid GN model and viscous NS model is good in most cases. Thus, this demonstrates that the effect of the viscosity on the generated waves is small. The oscillatory disturbance can produce periodic waves over the tank and, thus, can be designed as an innovative wave-maker or used in wave energy systems.