Internal solitary waves generated by a moving bottom disturbance

Abstract

The strongly nonlinear Miyata–Choi–Camassa model under the rigid lid approximation (MCC-RL model) can describe accurately the dynamics of large-amplitude internal waves in a two-layer fluid system for shallow configurations. In this paper, we apply the MCC-RL model to study the internal waves generated by a moving body on the bottom. For the case of the moving body speed $U=1.1c_{0}$, where ${c_0}$ is the linear long-wave speed, the accuracy of the MCC-RL results is assessed by comparing with Euler's solutions, and very good agreement is observed. It is found that when the moving body speed increases from $U=0.8c_{0}$ to $U=1.241c_{0}$, the amplitudes of the generated internal solitary waves in front of the moving body become larger. However, a critical moving body speed is found between $U=1.241c_{0}$ and $U=1.242c_{0}$. After exceeding this critical speed, only one internal wave right above the body is generated. When the moving body speed increases from $U=1.242c_{0}$ to $U=1.5c_{0}$, the amplitudes of the internal waves become smaller.

1. Introduction

Internal solitary waves exist widely in various seas all over the world; see e.g. Jackson (2007). Internal solitary waves are seen frequently in the northern part of the South China Sea, where Huang et al. (2016) observed an internal solitary wave with an amplitude as large as $240$ m.

To describe the internal solitary waves, the two-layer fluid system with constant mass densities is always considered. The Korteweg–de Vries (KdV) model is the earliest and widely used approach because of its simple form; see e.g. Benjamin (1966) and Miles (1980). However, when the KdV model describes large-amplitude internal solitary waves, the wave profiles are narrower and the wave speed is obviously faster when compared to the experimental data (Grue et al. 1999; Kodaira et al. 2016). Thus the KdV model is not suitable for describing large-amplitude internal waves due to the assumption of weak nonlinearity (Ostrovsky & Stepanyants 2005; Helfrich & Melville 2006).

Considering the large-amplitude internal waves in a two-layer shallow configuration (i.e. $h_1/\lambda \ll 1$ and $h_2/\lambda \ll 1$, where $h_1$ and $h_2$ are the depths of the upper and lower fluid layers, respectively, and $\lambda$ is the characteristic wavelength), Miyata (1985, 1988) and Choi & Camassa (1999) (MCC) derived a strongly nonlinear internal wave model. In their model, the depth-averaged horizontal velocities are used to describe the horizontal velocity variations along the fluid column for the upper and lower fluid layers. Here, we refer to it as the MCC-RL model because the free surface is approximated as a rigid lid (RL). Choi & Camassa (1996) also derived the MCC model that included the free surface effect (MCC-FS model).

Because of the excellent performance in describing the large-amplitude internal waves, the MCC model is widely used in various internal wave problems. These problems can be divided into three types based on the condition of the bottom boundary, namely (i) the flat bottom ($z=-h_{2}$), (ii) the space-varying bottom ($z=-h_{2}+b(x)$, where $b$ is the elevation of the bottom), and (iii) the time-varying bottom ($z=-h_{2}+b(x,t)$). Table 1 shows some related research on the applications of the MCC model to these three types of problems.

Table 1. Available literature on the application of the MCC model for the three types of bottom boundary conditions discussed in the text.

For the flat bottom problems, the MCC model is used widely to describe the large-amplitude internal waves in a two-layer fluid system (Choi & Camassa 1999; Camassa et al. 2006; Barros & Gavrilyuk 2007; Jo & Choi 2008; Kodaira et al. 2016; la Forgia & Sciortino 2019). Recently, Barros et al. (2020) extended the two-layer MCC-RL model to the three-layer case to study the properties of large-amplitude mode-2 internal solitary waves. Also, la Forgia & Sciortino (2020, 2021) used the MCC-RL model and the MCC-FS model, respectively, to study the internal solitary waves in the presence of a uniform current. Choi (2022) derived a second-order model to include the next-order correction of the MCC model. Good agreement was found between the results provided by the second-order model derived by Choi (2022) and the Euler's solutions.

For the space-varying bottom problems, Jo & Choi (2002) studied the deformation of an elevation internal solitary wave propagating over topography by use of the MCC-RL model, and the MCC-RL results were shown to be different from the weakly nonlinear prediction. Choi et al. (2020) used the MCC-RL model to study the propagation of a depression internal solitary wave over an isolated bottom topography, and a Fourier filter was applied to eliminate the local instability.

For the time-varying bottom problems, the MCC model has been extended for the case of a time-varying bottom (and also multiple layers) by Choi (2000). It is expected that the MCC model can be applied to provide further information about the time-varying bottom problems. Internal waves generated by a moving disturbance have been a subject of great interest. Internal waves may be generated by a surface disturbance, such as the dead-water phenomenon induced by ship motion (see e.g. Mercier, Vasseur & Dauxois 2011; Duchene 2011) or by a seafloor disturbance, such as underwater landslides (see e.g. Brizuela, Filonov & Alford 2019). In the present study, we focus on internal waves generated by a moving body on the seafloor. Under the rigid lid assumption, Grue et al. (1997) established a time-stepping method for solving Euler's equations to study this problem in a two-layer fluid system. They showed that a moving body with speed $U=1.1c_{0}$ (where ${c_0}$ is the linear long-wave speed) could generate a series of internal solitary waves. The number of waves would increase with the moving distance. However, the effect of the moving body speed on the generated internal waves should be investigated further.

The motivations of this study are (i) to apply the MCC-RL model to the time-varying bottom problems, and (ii) to analyse the effect of the moving body speed on the generated internal waves.

This paper is organized as follows. In § 2, the equations of the MCC-RL model with time-varying bottom are derived. In § 3, the numerical algorithm is presented. Numerical test cases are presented and discussed in § 4. Conclusions are reached in § 5.

2. The MCC-RL model with time-varying bottom

In this section, the MCC-RL model with time-varying bottom is introduced. We consider a two-dimensional system of two fluid layers whose densities and undisturbed thicknesses are given by $\rho _i$ and $h_i$, respectively, where $i=1$ represents the upper fluid layer and $i=2$ represents the lower fluid layer. The origin of the two-dimensional Cartesian coordinate system is set at the undisturbed interface between the two fluid layers, $x$ is the horizontal axis, positive to the right, and $z$ is the vertical axis, positive up. The upper surface of the upper fluid layer, the interface between the two fluid layers and the lower surface of the lower fluid layer are represented by $z=h_1$, $z=\zeta (x,t)$ and $z=-h_2+b(x,t)$, respectively, where $b(x,t)$ is the bottom elevation. Also, $\eta _1(x,t)$ and $\eta _2(x,t)$ are the local layer thicknesses of the upper fluid layer and lower fluid layer, respectively. A sketch of a two-layer fluid system where the bottom varies with time is shown in figure 1.

Figure 1. Sketch of a two-layer fluid system where the bottom varies with time.

For incompressible, homogeneous fluids, the mass conservation equations of the two fluid layers can be written as

(2.1a)$$\begin{gather} {u_{1,x}} + {w_{1,z}} = 0, \end{gather}$$

(2.1b)$$\begin{gather}{u_{2,x}} + {w_{2,z}} = 0, \end{gather}$$

where $u_{i}(x,z,t)$ and $w_{i}(x,z,t)$ are the horizontal velocity and vertical velocity respectively, where $i=1$ represents the upper fluid layer and $i=2$ represents the lower fluid layer. The subscripts $x$ and $z$ after a comma represent the spatial partial derivatives.

For this fluid system, and assuming that viscous effects are negligible, the momentum conservation equations for the two fluid layers can be written as

(2.2a)$$\begin{gather} {u_{1,t}} + {u_1}{u_{1,x}} + {w_1}{u_{1,z}} = {{ - {p_{1,x}}} / {{\rho _1}}}, \end{gather}$$

(2.2b)$$\begin{gather}{u_{2,t}} + {u_2}{u_{2,x}} + {w_2}{u_{2,z}} = {{ - {p_{2,x}}} / {{\rho _2}}}, \end{gather}$$

(2.2c)$$\begin{gather}{w_{1,t}} + {u_1}{w_{1,x}} + {w_1}{w_{1,z}} = {{ - {p_{1,z}}} / {{\rho _1}}} - g, \end{gather}$$

(2.2d)$$\begin{gather}{w_{2,t}} + {u_2}{w_{2,x}} + {w_2}{w_{2,z}} = {{ - {p_{2,z}}} / {{\rho _2}}} - g, \end{gather}$$

where $p_{1}(x,z,t)$ and $p_{2}(x,z,t)$ are the pressures of the upper fluid layer and the lower fluid layer, respectively, and $g$ is the gravitational acceleration. The subscript $t$ after a comma indicates the partial derivative with respect to time.

The kinematic boundary conditions for the upper fluid layer are written as

(2.3a)$$\begin{gather} {w_1} = 0,\quad{\rm at}\ z = {h_1}, \end{gather}$$

(2.3b)$$\begin{gather}{w_1} = {\zeta _{,t}} + {u_1}{\zeta _{,x}},\quad {\rm at}\ z = {\zeta}(x,t). \end{gather}$$

The kinematic boundary conditions for the lower fluid layer are written as

(2.4a)$$\begin{gather} {w_2} = {\zeta_{,t}} + {u_2}{\zeta_{,x}},\quad {\rm at}\ z = {\zeta}(x,t), \end{gather}$$

(2.4b)$$\begin{gather}{w_2} = {b_{,t}} + {u_2}{b_{,x}}, \quad {\rm at}\ z =-h_2+b(x,t). \end{gather}$$

The dynamic boundary condition at the interface between the two fluid layers is written as

(2.5)\begin{equation} {\mathop{p}\limits^{\smallsmile}}_{1} = {\mathop{p}\limits^{\smallfrown}}_{2}=P,\quad {\rm at}\ z = {\zeta}(x,t), \end{equation}

where ${{\mathop {p}\limits ^{\smallsmile }}_{1}(x,t)}$ is the pressure at the lower surface of the upper fluid layer, ${{\mathop {p}\limits ^{\smallfrown }}_{2}(x,t)}$ is the pressure at the upper surface of the lower fluid layer and $P(x,t)$ is the pressure at the interface.

In the MCC-RL model, it is assumed that the characteristic wavelength is long compared with each fluid layer thickness. When the bottom is flat, i.e. $b(x,t)=0$, the equations of the MCC-RL model given by Choi & Camassa (1999) can be written in terms of four unknowns ($\zeta,{\bar u_1},{\bar u_{2}},P$) as

(2.6a)$$\begin{gather} {\eta _{1,t}} + {\left( {{\eta _1}{{\bar u}_1}} \right)_{,x}} = 0, \quad{\eta _1} = {h_1} - \zeta, \end{gather}$$

(2.6b)$$\begin{gather}{\eta _{2,t}} + {\left( {{\eta _2}{{\bar u}_2}} \right)_{,x}} = 0, \quad {\eta _2} = {h_2} + \zeta, \end{gather}$$

(2.6c)$$\begin{gather}{{\bar u}_{1,t}} + {{\bar u}_1}{{\bar u}_{1,x}} + g{\zeta _{,x}} =- \frac{{{P_{,x}}}}{{{\rho _1}}} + \frac{1}{{{\eta _1}}}{\left( {\frac{1}{3}\,\eta_1^3{G_1}} \right)_{,x}}, \end{gather}$$

(2.6d)$$\begin{gather}{{\bar u}_{2,t}} + {{\bar u}_2}{{\bar u}_{2,x}} + g{\zeta _{,x}} =- \frac{{{P_{,x}}}}{{{\rho _2}}} + \frac{1}{{{\eta _2}}}{\left( {\frac{1}{3}\,\eta_2^3{G_2}} \right)_{,x}}, \end{gather}$$

where ${\bar u_1}$ and ${\bar u_2}$ are the depth-averaged horizontal velocities, which are defined as

(2.7a)$$\begin{gather} {\bar u_1}( {x,t} ) = \frac{1}{{{\eta _1}}}\int_\zeta ^{{h_1}} {{u_1}( {x,z,t} )\,\text{d}z}, \end{gather}$$

(2.7b)$$\begin{gather}{\bar u_2}( {x,t} ) = \frac{1}{{{\eta _2}}}\int_{ - {h_2}}^\zeta {{u_2}( {x,z,t} )\,\text{d}z}, \end{gather}$$

and $G_{1}$ and $G_{2}$ are defined as

(2.8a,b)\begin{equation} {G_1}( {x,t} ) = {\bar u_{1,xt}} + {\bar u_1}{\bar u_{1,xx}} - {\left( {{{\bar u}_{1,x}}} \right)^2} ,\quad {G_2}( {x,t} ) = {\bar u_{2,xt}} + {\bar u_2}{\bar u_{2,xx}} - {\left( {{{\bar u}_{2,x}}} \right)^2}. \end{equation}

When the bottom varies with time and space, by the multi-layer MCC model proposed by Choi (2000), we obtain the equations for the MCC-RL model for the four unknowns ($\zeta,{\bar u_1},{\bar u_2},P$) as

(2.9a)$$\begin{gather} {\eta _{1,t}} + {\left( {{\eta _1}{{\bar u}_1}} \right)_{,x}} = 0, \quad {\eta _1} = {h_1} - \zeta, \end{gather}$$

(2.9b)$$\begin{gather}{\eta _{2,t}} + {\left( {{\eta _2}{{\bar u}_2}} \right)_{,x}} = 0, \quad {\eta _2} = {h_2} + \zeta - b, \end{gather}$$

(2.9c)$$\begin{gather}{{\bar u}_{1,t}} + {{\bar u}_1}{{\bar u}_{1,x}} + g{\zeta _{,x}} =- \frac{{{P_{,x}}}}{{{\rho _1}}} + \frac{1}{{{\eta _1}}}{\left( {\frac{1}{3}\,\eta_1^3{G_1}} \right)_{,x}}, \end{gather}$$

(2.9d)$$\begin{gather} {{\bar u}_{2,t}} + {{\bar u}_2}{{\bar u}_{2,x}} + g{\zeta _{,x}} =- \frac{{{P_{,x}}}}{{{\rho _2}}} + \frac{1}{{{\eta _2}}}{\left( {\frac{1}{3}\,\eta_2^3{G_2} - \frac{1}{2}\,\eta_2^2D_2^2b} \right)_{,x}}\nonumber\\+ \left( {\frac{1}{2}\,{\eta _2}{G_2} - D_2^2b} \right){b_{,x}}, \end{gather}$$

where ${D_2} \equiv \partial t + {\bar u_2}\,\partial x$. For a given bottom elevation $b$, (2.9) are closed and are solvable.

Considering that the bottom varies with time and space, there are some differences between (2.9) and (2.6). In (2.9b), compared with (2.6b), the local layer thickness $\eta _2$ is changed. In (2.9d), compared with (2.6d), some terms related to time-varying bottom are added, including $- ({1}/{{2{\eta _2}}}){( {\eta _2^2D_2^2b} )_{,x}}$ and $( {\tfrac {1}{2}{\eta _2}{G_2} - D_2^2b} ){b_{,x}}$.

In addition, considering the bottom elevation, compared with (2.7b), ${\bar u_2}$ is defined as

(2.10)\begin{equation} {\bar u_2}( {x,t} ) = \frac{1}{{{\eta _2}}}\int_{ - {h_2}+b}^\zeta {{u_2}( {x,z,t} )\,\text{d}z}. \end{equation}

By eliminating $\zeta _{,t}$ from (2.9a) and (2.9b), ${\bar u_1}$ can be expressed in terms of ${\bar u_2}$, $\zeta$ and $b$ as

(2.11)\begin{equation} {\bar u_1} = \frac{1}{{{\eta _1}}}\left( {\int_{ - \infty }^{ + \infty } {{b_{,t}}\,\text{d} x} - {{\bar u}_2}{\eta _2}} \right). \end{equation}

3. Numerical algorithm

In the time domain algorithm for solving the MCC-RL model, Jo & Choi (2008) used the second-order central difference scheme both in space and in time. In this section, we will give a numerical algorithm with higher accuracy for solving the MCC-RL equations.

For the equations of the MCC-RL model with time-varying bottom, combining (2.9c) and (2.9d) to eliminate $P_{,x}$, we can obtain

(3.1)$$\begin{gather} \frac{{{\rho _1}}}{{{\eta _1}}}{\left( {\frac{1}{3}\,\eta_1^3{G_1}} \right)_{,x}} - {\rho _1}\left( {{{\bar u}_{1,t}} + {{\bar u}_1}{{\bar u}_{1,x}} + g{\zeta _{,x}}} \right) = \frac{{{\rho _2}}}{{{\eta _2}}}{\left( {\frac{1}{3}\,\eta_2^3{G_2} - \frac{1}{2}\,\eta_2^2D_2^2b} \right)_{,x}}\nonumber\\ {}+ {\rho _2}\left( {\frac{1}{2}\,{\eta _2}{G_2} - D_2^2b} \right){b_{,x}} - {\rho _2}\left( {{{\bar u}_{2,t}} + {{\bar u}_2}{{\bar u}_{2,x}} + g{\zeta _{,x}}} \right), \end{gather}$$

where ${\bar u_1}$ can be eliminated by using (2.11).

Then (3.1) can be arranged in the form

(3.2)\begin{equation} A{\bar u_{2,xxt}} + B{\bar u_{2,xt}} + C{\bar u_{2,t}} = F, \end{equation}

where $A$, $B$, $C$ and $F$ are functions of $\zeta (x,t)$, $b(x,t)$, $\bar u_{2}(x,t)$ and their spatial derivatives. Here, $A$, $B$, $C$ and $F$ are defined as

(3.3a)\begin{align} A &=-\frac{1}{3}\,{\eta _2}\left( {{\rho _1}{\eta _1} + {\rho _2}{\eta _2}} \right), \end{align}

(3.3b)\begin{align} B& =- \frac{1}{3}\left[ {{\rho _1}{\eta _2}{\eta _{1,x}} + \left( {2{\rho _1}{\eta _1} + 3{\rho _2}{\eta _2}} \right){\eta _{2,x}}} \right], \end{align}

(3.3c)\begin{align} C &= \frac{{{\rho _1}}}{3}\left( {{\eta _{1,xx}}{\eta _2} - {\eta _{1,x}}{\eta _{2,x}} - \eta _1^{}{\eta _{2,xx}} + {\eta _2}\,\frac{{3 + \eta _{1,x}^2}}{{{\eta _1}}}} \right)\nonumber\\ &\quad + {\rho _2}\left( {1 + b_{,x}^2 + {b_{,x}}{\eta _{2,x}} + \frac{{{b_{,xx}}{\eta _2}}}{2}} \right), \end{align}

(3.3d)\begin{align} F& = A{{\bar u}_{2,xxt}} + B{{\bar u}_{2,xt}} + C{{\bar u}_{2,t}} - \frac{{{\rho _1}}}{{{\eta _1}}}{\left( {\frac{1}{3}\,\eta_1^3{G_1}} \right)_{,x}}\nonumber\\ & \quad + {\rho _1}\left( {{{\bar u}_{1,t}} + {{\bar u}_1}{{\bar u}_{1,x}} + g{\zeta _{,x}}} \right) + \frac{{{\rho _2}}}{{{\eta _2}}}{\left( {\frac{1}{3}\,\eta_2^3{G_2} - \frac{1}{2}\,\eta_2^2D_2^2b} \right)_{,x}}\nonumber\\ &\quad + {\rho _2}\left( {\frac{1}{2}\,{\eta _2}{G_2} - D_2^2b} \right){b_{,x}} - {\rho _2}\left( {{{\bar u}_{2,t}} + {{\bar u}_2}{{\bar u}_{2,x}} + g{\zeta _{,x}}} \right). \end{align}

The spatial difference discretization is utilized to solve (3.2). The calculation domain of $x$ is discretized into having a uniform grid of $x$ values, spaced $\Delta {x}$ apart. The $i$th point on the grid is denoted by ${x_i} = i\,\Delta {x}$ for $i = 1, 2, \ldots, n$. Time is discretized with intervals $\Delta {t}$, with ${t_j} = j\,\Delta {t}$. For example, the value $\bar u_{2}(x_i,t_j)$ will be denoted by $\bar u_2^{( i )}$, where $j$ is omitted since we refer to the same $j$ time. Similar subscripts are also used for other variables.

The five-point central difference scheme is used to calculate spatial derivations of $\bar u_{2,xt}^{( {i } )}$ and $\bar u_{2,xxt}^{( {i } )}$ as follows:

(3.4a)$$\begin{gather} \bar u_{2,xt}^{(i)} = \frac{{\bar u_{2,t}^{(i-2)} - 8\bar u_{2,t}^{(i-1)} + 8\bar u_{2,t}^{(i+1)} - \bar u_{2,t}^{(i+2)}}}{{12\Delta x}}, \end{gather}$$

(3.4b)$$\begin{gather}\bar u_{2,xxt}^{(i)} = \frac{{ - \bar u_{2,t}^{(i-2)} + 16\bar u_{2,t}^{(i-1)} - 30\bar u_{2,t}^{(i)} + 16\bar u_{2,t}^{(i+1)} - \bar u_{2,t}^{(i+2)}}}{{12{{\left( {\Delta x} \right)}^2}}}. \end{gather}$$

Substituting (3.4) into (3.2) will result in

(3.5)\begin{equation} {\tilde A^{(i)}}\bar u_{2,t}^{(i-2)} + {\tilde B^{(i)}}\bar u_{2,t}^{(i-1)} + {\tilde C^{(i)}}\bar u_{2,t}^{(i)} + {\tilde D^{(i)}}\bar u_{2,t}^{(i+1)} + {\tilde E^{(i)}}\bar u_{2,t}^{(i+2)} = {F^{(i)}}, \end{equation}

where

(3.6a)$$\begin{gather} {\tilde A^{(i)}} =- {A^{(i)}}\,\frac{1}{{12{{\left( {\Delta x} \right)}^2}}} + {B^{(i)}}\,\frac{1}{{12\Delta x}}, \end{gather}$$

(3.6b)$$\begin{gather}{\tilde B^{(i)}} = { A ^{(i)}}\,\frac{{4}}{{3{{\left( {\Delta x} \right)}^2}}} - { B ^{(i)}}\,\frac{2}{{3\Delta x}}, \end{gather}$$

(3.6c)$$\begin{gather}{\tilde C^{(i)}} =- {A ^{(i)}}\,\frac{{5}}{{2{{\left( {\Delta x} \right)}^2}}} + {C ^{(i)}}, \end{gather}$$

(3.6d)$$\begin{gather}{\tilde D^{(i)}} = {A ^{(i)}}\,\frac{{4}}{{3{{\left( {\Delta x} \right)}^2}}} + { B ^{(i)}}\,\frac{2}{{3\Delta x}}, \end{gather}$$

(3.6e)$$\begin{gather}{\tilde E^{(i)}} =- { A ^{(i)}}\,\frac{1}{{12{{\left( {\Delta x} \right)}^2}}} -{ B ^{(i)}}\,\frac{1}{{12\Delta x}}. \end{gather}$$

Details about the algorithm used to solve (3.5) to obtain $\bar u_{2,t}$ can be found in Zhao, Duan & Ertekin (2014). Meanwhile, $\zeta _{,t}$ can be calculated by (2.9b). We use the fourth-order Adams predictor–corrector scheme for time marching. For each case, we have tested the convergence of $\Delta x$ and $\Delta t$. Here, the results shown are the converged results.

In the time domain simulation of internal waves in a two-layer system, the time-dependent inviscid model suffers from the Kelvin–Helmholtz instability due to the velocity discontinuity across the interface, as mentioned by Jo & Choi (2002). In order to reduce the effect of local instability, a numerical filter is found to be effective to suppress the short-wave instability without affecting the long-wavelength behaviour; see e.g. Jo & Choi (2008) and Choi et al. (2020). On the other hand, the internal wave model can be regularized to eliminate shear instability; see e.g. Choi, Barros & Jo (2009), Lannes & Ming (2015) and Duchêne, Israwi & Talhouk (2016).

In our simulation, a five-point smoothing filter is applied intermittently in the numerical solutions for the time-stepping variables $\zeta$ and $\bar u_{2}$ to dampen the high wavenumber disturbances, and this filter is written as

(3.7a)$$\begin{gather} {\dot f^{(i)}} = \frac{{17{f^{(i)}} + 12\left( {{f^{(i-1)}} + {f^{(i+1)}}} \right) - 3\left( {{f^{(i-2)}} + {f^{(i+2)}}} \right)}}{{35}}, \end{gather}$$

(3.7b)$$\begin{gather}{\ddot f^{(i)}} = {\dot f^{(i)}}*\gamma + {f^{(i)}}*\left( {1 - \gamma } \right), \end{gather}$$

where ${{\dot f}^{( i )}}$ is the variable after the first step of smoothing, ${{\ddot f}^{( i )}}$ is the final variable after smoothing and $\gamma$ is the smoothing parameter used for weighting process, typically taken as 0.01. The smoothing formula has a minor effect on numerical results (Longuet-Higgins & Cokelet 1976; Fuhrman, Madsen & Bingham 2006; Zhao et al. 2014).

In order to test the effect of the smoothing filter, we simulate the internal solitary wave that suffers from Kelvin–Helmholtz instability, shown in Jo & Choi (2002). The parameters are $\rho _1/\rho _2=1/1.01$, $h_1/h_2=1/2$ and $a/h_1=-0.4885$. As shown in figure 2, the smoothing filter in our simulations can effectively reduce the influence of the short waves and ensure stable propagation of large-amplitude internal waves over a long time. By comparing the wave profiles at $t(g/h_{1})=0$ with $t(g/h_{1})=4000$, we find that the results are almost identical and there is no numerical dissipation.

Figure 2. The propagation of a large-amplitude internal solitary wave over an extended time in a frame moving with the internal solitary wave speed, with $\rho _1/\rho _2=1/1.01$, $h_1/h_2=1/2$ and $a/h_1=-0.4885$. (a) The propagation of the internal solitary wave in the space–time domain. (b) Internal solitary wave profiles.

4. Numerical test cases

In this section, we will use the MCC-RL model to conduct a number of numerical tests. In the first case, we consider a large-amplitude internal solitary wave propagating on the flat bottom and compare the wave profile and velocity field obtained by the model discussed here with the experimental data and Euler's solutions given by Grue et al. (1999). In the second case, we study the internal waves generated by the moving body with speed $U=1.1c_{0}$ on the bottom, and compare the MCC-RL results with Euler's solutions given by Grue et al. (1997). Next, we study the effect of moving body speed on the generated internal waves, with the speed ranging from $U=0.8c_{0}$ to $U=1.5c_{0}$. Finally, we apply the MCC-RL model to simulate the internal waves generated by an unsteady moving bottom.

4.1. Internal solitary wave propagating on the flat bottom

In this subsection, we consider a large-amplitude internal solitary wave propagating on the flat bottom. Following the physical experiments conducted by Grue et al. (1999), we select the parameters as follows: $h_1 =0.15$ m, $h_2 =0.62$ m, $\rho _1 =999\,\text {kg}\,\text {m}^{-3}$ and $\rho _2 =1022\,\text {kg}\,\text {m}^{-3}$. The amplitude of the internal solitary wave that we selected is $-1.23{h_1}$. The initial values are provided by the steady solutions of the MCC-RL model. More details on the steady solution of MCC-RL model can be found in Choi & Camassa (1999).

A snapshot of the internal solitary wave propagating over the flat bottom at different moments is shown in figure 3(a). We translate the internal wave profiles at different moments to the place where the crest is at $x/h_1=0$ in figure 3(b). From $t=0$ s to $t=600$ s, we find that the profiles of the propagating internal wave show very good agreement, which indicates that this large-amplitude internal solitary wave can propagate on the flat bottom steadily based on the numerical algorithm discussed in § 3.

Figure 3. A large-amplitude internal solitary wave propagating on the flat bottom in a frame moving with the internal solitary wave speed, with $\rho _1/\rho _2=0.977$, $h_1/h_2=1/4.13$ and $a/h_1=-1.23$. (a) The propagation of the internal solitary wave. (b) Translated internal solitary wave profiles at different moments (the lines are almost exactly on top of each other).

In figure 4, we present directly the time domain results on the internal wave profile and velocity field at $t=600$ s. Euler's solutions and experimental data obtained by Grue et al. (1999) are also shown for comparison purposes. For the horizontal velocity along the fluid column at the internal wave crest in figure 4(b), the MCC-RL result can be obtained based on the depth-averaged horizontal velocities ${\bar u_1}$ and ${\bar u_2}$ (Camassa et al. 2006). In figure 4(b), $c_0$ is the linear long-wave speed, which is defined as

(4.1)\begin{equation} {c_0} = \sqrt {{{g{h_1}{h_2}\left( {{\rho _2} - {\rho _1}} \right)}/{\left( {{\rho _1}{h_2} + {\rho _2}{h_1}} \right)}}}. \end{equation}

Figure 4. Comparison between the time domain solution of the MCC-RL model and Euler's solution and experimental data, with $\rho _1/\rho _2=0.977$, $h_1/h_2=1/4.13$ and $a/h_1=-1.23$. (a) Internal wave profile. (b) Horizontal velocity along the fluid column at the internal wave crest.

From figure 4, we find that the numerical results show good agreement with each other, and they both match the experimental data very well. Thus it is demonstrated that the time domain solution of the MCC-RL model in this case is accurate.

4.2. Internal waves generated by a moving bottom disturbance

In this subsection, we consider the internal waves generated by a moving body on the bottom. The continuous and uniform motion of a body on the bottom, in a moving coordinate system, could also be regarded as a tidal flow over topography due to the relative motions (Wang 2019). Various approaches are used in the literature to investigate this problem, including by use of the forced KdV model in a single-layer fluid (Grimshaw 2010; Grimshaw & Maleewong 2016) and in a continuously density-stratified fluid (Grimshaw & Smyth 1986; Grimshaw, Chan & Chow 2002; Grimshaw & Helfrich 2018). As the forced KdV model is a weakly nonlinear model with limitations, it is desirable to use a strongly nonlinear model to further study the internal waves generated by a moving body on the bottom (Helfrich & Melville 2006).

In the simulation, we use the same parameters given by Grue et al. (1997), as follows: $h_1 =0.12$ m, $h_2 =0.03$ m, $\rho _1 =787.3\,\text {kg}\,\text {m}^{-3}$ and $\rho _2 =1000\,\text {kg}\,\text {m}^{-3}$. The moving bottom disturbance is a semi-ellipse (shown in figure 5), whose shape is fixed at all times, and it is confined to only horizontal motion on the seafloor. The semi-major axis and the semi-minor axis of the moving ellipse body are ${L_{1/2}} = 10{h_2} = 0.3$ m and ${B_{1/2}} = 0.1{h_2} = 0.003$ m, respectively. The constant speed of the body is $U = 1.1{c_0} = 0.252\,{\rm m}\,{\rm s}^{-1}$, where the value of $c_{0}$ is obtained by (4.1) as $0.229\,{\rm m}\,{\rm s}^{-1}$. It should be noted that in this study, the semi-ellipse on the bottom is moving with the constant speed starting from the initial moment, i.e. there is no acceleration of the bottom disturbance. This, however, is not required in general. A sketch of the physical problem is shown in figure 5.

Figure 5. Sketch of semi-ellipse moving on the bottom in a two-layer fluid system.

The internal surface elevations at different moments obtained by the MCC-RL model are shown in figure 6, where Euler's solutions obtained by Grue et al. (1997) are also presented for comparison.

Figure 6. Internal elevations generated by a moving body on the bottom at different moments, with $\rho _1/\rho _2=0.7873$, $h_1/h_2=4/1$ and $U=1.1c_{0}$, for (a) $t=59.7$ s, (b) $t=106.2$ s, (c) $t=152.7$ s.

As shown in figure 6, several internal waves can be generated in front of the moving body, and with the moving time increasing, the number of generated internal waves increases. By comparing the MCC-RL results with the Euler's solutions obtained by Grue et al. (1997), we observe that the two results show good agreement in general since this case belongs to the shallow configuration case. We also observe in figure 6(c) that the amplitudes of the generated internal waves obtained by the MCC-RL model are approximately $0.84h_{2}$.

In figure 7, we compare the time domain solution at $152.7$ s with the steady solution of the MCC-RL model on internal wave profiles. The steady solution of the MCC-RL model can be obtained by (3.50) and (3.57) of Choi & Camassa (1999), as are also shown below:

(4.2a)$$\begin{gather} {\left( {{\zeta _{,X}}} \right)^2} = \frac{{3{\zeta ^2}\left[ {{\rho _1}c_{{w}}^2{\eta _2} + {\rho _2}c_{{w}}^2{\eta _1} - g\left( {{\rho _2} - {\rho _1}} \right){\eta _1}{\eta _2}} \right]}}{{{\rho _1}c_{{w}}^2h_1^2{\eta _2} + {\rho _2}c_{{w}}^2h_2^2{\eta _1}}}, \end{gather}$$

(4.2b)$$\begin{gather}\frac{{c_{{w}}^2}}{{c_0^2}} = \frac{{\left( {{h_1} - a} \right)\left( {{h_2} + a} \right)}}{{{h_1}{h_2} - \left( {{{c_0^2} / g}} \right)a}}, \end{gather}$$

where $X = x - c_{{w}}t$, $c_{{w}}$ is the speed of the internal solitary wave and $a$ is the amplitude of the internal solitary wave. Good agreements are observed, and this indicates that the generated waves are indeed internal solitary waves. Furthermore, it is found that the distances between the crests of adjacent internal solitary waves are approximately $27h_{2}$. This indicates that there is a weak interaction between the internal solitary waves. Similar conclusions are also made by Grue et al. (1997).

Figure 7. Comparison of the internal wave profiles between the time domain solution and the steady solution by the MCC-RL model, with $\rho _1/\rho _2=0.7873$, $h_1/h_2=4/1$ and $a/h_2=0.84$.

In table 2, the positions of the first three wave crests at $t=150$ s and $180$ s are given. According to table 2, the speeds of these waves is approximately $1.23c_{0}$, which is the same as the speed provided by the steady solution. Thus the generated internal waves can be regarded as a series of internal solitary waves.

Table 2. Positions of the wave crests.

4.3. Effect of the moving body speed on the generated internal waves

In this subsection, we focus on the effect of the moving body speed on the generated internal waves. By conducting a series of numerical simulations, we find that when the moving body speed is $U<0.8c_{0}$, the amplitudes of generated internal waves are quite small. Thus the constant moving body speeds that we select are from $U=0.8{c_0}$ to $U=1.5{c_0}$, namely $U=0.8{c_0}$, $1.0{c_0}$, $1.241{c_0}$, $1.242{c_0}$, $1.4{c_0}$ and $1.5{c_0}$. The initial centre position of the moving body is located at $x_{0}=5$ m. Other parameters are the same as those given in § 4.2. For convenience of display, we translate the moving body centre to $x=0$ in the following figures.

4.3.1. $U=0.8{c_0}$

For the case of the moving body speed $U=0.8{c_0}$, we compare the generated internal waves at $t=60$ s, $120$ s and $180$ s, shown in figure 8(a). Comparing with the case of the moving body speed $U=1.1{c_0}$ shown in figure 6, we find that the amplitude of the generated internal waves is smaller in this case. The amplitude of the leading internal wave is only $0.22h_{2}$. In figure 8(b), we compare the time domain solution at $t=180$ s with the steady solution on the first two internal wave profiles, and good agreements are observed. Moreover, the amplitude of the generated second internal wave is smaller than that of the leading internal wave due to their interaction.

Figure 8. Generation of internal waves with the moving body speed $U=0.8c_{0}$, with $\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$. (a) Comparison of the internal wave profiles at different moments. (b) Comparison between the time domain solution and the steady solution.

The amplitudes of the first two internal waves and the crest distance at different moments are shown in table 3. From $t=60$ s to $t=180$ s, the amplitude of the leading wave $a_{1}$ is basically unchanged, while the amplitude of the second wave $a_{2}$ increases. Meanwhile, the distance between the wave crests, $l_{1-2}$, also becomes larger when the time increases, which indicates that their interaction becomes weaker.

Table 3. The amplitudes and positions of the first two generated internal waves at different moments for the case $U=0.8c_{0}$, where ${a_1}$ is the amplitude of the leading wave, ${a_2}$ is the amplitude of the second wave and ${l_{1 - 2}}$ is the crest distance between the leading wave and the second wave.

4.3.2. $U=1.0{c_0}$

For the case of the moving body speed $U=1.0{c_0}$, the generated internal waves at $t=60$ s, $120$ s and $180$ s are shown in figure 9(a), where we observe that the amplitudes of the generated internal wave become larger than those of the case of the moving body speed $U=0.8{c_0}$. The amplitude of the leading generated wave reaches $0.58h_{2}$. Good agreements are found between the time domain solution at $t=180$ s and the steady solution for the first two generated internal wave profiles in general.

Figure 9. Generation of internal waves with the moving body speed $U=1.0c_{0}$, with $\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$. (a) Comparison of the internal wave profiles at different moments. (b) Comparison between time domain solution and steady solution.

The amplitudes of the first two generated internal waves and the crest distance at different moments are shown in table 4. From $t=60$ s to $t=180$ s, the amplitudes of the generated leading internal wave are basically unchanged, and the amplitudes of the second wave become slightly larger. Also, the distance between the wave crests $l_{1-2}$ does not show obvious differences when the time increases, which indicates that there is weak interaction between the two internal waves.

Table 4. The amplitudes and positions of the first two generated internal waves at different moments for the case $U=1.0c_{0}$, where ${a_1}$ is the amplitude of the leading wave, ${a_2}$ is the amplitude of the second wave and ${l_{1 - 2}}$ is the crest distance between the leading wave and the second wave.

4.3.3. $U=1.241{c_0}$

For the case of the moving body speed $U=1.241{c_0}$, the amplitudes of the generated internal waves increases further, and the time required to generate the internal wave needs longer. Thus we compare the generated internal waves at $t=140$ s, $280$ s and $420$ s in figure 10(a). At $t=420$ s, the amplitude of the leading internal wave reaches $1.23h_{2}$. In figure 10(b), good agreement can be found between the time domain results and the steady solution.

Figure 10. Generation of internal waves with the moving body speed $U=1.241c_{0}$, with $\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$. (a) Comparison of the internal wave profiles at different moments. (b) Comparison between time domain solution and steady solution.

Next, we focus on the relationship between the amplitude of the generated leading internal solitary wave and the moving body speed. We recall that the moving body speed varies between $U=0.8c_{0}$ and $U=1.241c_{0}$. As shown in figure 11, the amplitude of the leading internal wave increases monotonically with an increase in the moving body speed. A similar phenomenon was also observed in the experiments conducted by Melville & Helfrich (1987), although the moving body was on the free surface in their experiments.

Figure 11. Relationship between leading internal wave amplitude $a_{1}$ and the moving body speed $U$, the moving speed changing from $U=0.8c_{0}$ to $U=1.241c_{0}$, with $\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$.

Furthermore, the relationship between the leading internal wave speed and the moving body speed is shown in figure 12. We find that the wave speed relative to the moving body speed, $c_{{w}}-U$, becomes smaller with the increasing moving body speed. When the moving body speed is $U=1.241{c_0}$, $c_{{w}}-U$ is only $0.037{c_0}$. This also explains that when the moving body speed is $U=1.241{c_0}$, it takes a longer time to generate internal waves.

Figure 12. Relationship between the generated internal wave relative speed $c_{{w}}-U$ and the moving body speed $U$, the moving speed changing from $U=0.8c_{0}$ to $U=1.241c_{0}$, with $\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$.

4.3.4. $U=1.242{c_0}$

For the case of the moving body speed $U=1.242{c_0}$ (and larger speeds), we find that the results are quite different from the results of the previous cases. For this case, we find that a single internal wave, right above the moving body, is generated. Apart from the time domain simulation, we also obtain the steady solution with time-varying bottom for comparison purposes.

Here, we introduce briefly the algorithm for obtaining the steady solution for the internal wave with a time-varying bottom. The speed of the moving reference frame is the same as the moving body speed $U$. We use the moving coordinates, which are located at the still interface of the two fluid layers, $XOZ$, to solve this problem. In the moving coordinates, $\zeta$, $b$, $\bar {u}_1$ and $\bar {u}_2$ can be written as

(4.3a)$$\begin{gather} \zeta (x,t) = \zeta (X), \end{gather}$$

(4.3b)$$\begin{gather}b(x,t) = b(X), \end{gather}$$

(4.3c)$$\begin{gather}{\bar u_1}(x,t) = {\bar u_1}(X), \end{gather}$$

(4.3d)$$\begin{gather}{\bar u_2}(x,t) = {\bar u_2}(X), \end{gather}$$

where $X = x - Ut$.

After substituting (4.3) into (2.9) and using the relation

(4.4)\begin{equation} {f_{,t}} =- {U}{f_{,X}}, \end{equation}

where $f=(\zeta, b, \bar u_1, \bar u_2)$, the moving coordinates form of the MCC-RL equations that are used to obtain the steady solution with time-varying bottom can be obtained. For a given bottom $b(X)$ and speed $U$, the system of equations is closed and solvable. The central difference scheme is used to calculate the spatial derivatives, and the Newton–Raphson method is used to determine the travelling solution. More details of the numerical scheme to solve the steady problem are given in Zhao et al. (2016) and Duan et al. (2018). For a given speed $U$, we can obtain the steady solution of the internal wave. It should be noted that it is only when the speed $U$ is large, such as $1.242c_{0}$, that the steady solution with time-varying bottom can be obtained. When the speed $U$ is less than the critical speed, a series of internal solitary waves will be generated continuously in front of the moving body, therefore steady state cannot be realized.

In figure 13, there is only one internal wave above the moving body, and the wave speed is the same as the moving body speed. In figure 13(a), the amplitude of the internal wave at $t=100$ s is slightly smaller than those at $t=200$ s and $300$ s. At $t=300$ s, we find that the amplitude of the generated internal wave is $0.41h_{2}$, which is obviously smaller than the $1.23h_2$ obtained in the case of the moving body speed $U=1.241c_0$. As shown in figure 13(b), good agreement can be found between the time domain solution at $t=300$ s and the steady solution with a time-varying bottom.

Figure 13. Generation of internal waves with the moving body speed $U=1.242c_{0}$, with $\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$. (a) Comparison of the internal wave profiles at different moments. (b) Comparison between time domain solution and steady solution.

4.3.5. $U=1.4{c_0}$

In figure 14, for the moving body speed $U=1.4{c_0}$, the amplitudes of the internal waves at $t=100$ s, $200$ s and $300$ s are basically the same. Comparing with the case of the moving body speed $U=1.242{c_0}$, the internal wave amplitude decreases to $0.17h_{2}$. As shown in figure 14(b), the time domain solution and the steady solution of the time-varying bottom are also in good agreement.

Figure 14. Generation of internal waves with the moving body speed $U=1.4c_{0}$, with $\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$. (a) Comparison of the internal wave profiles at different moments. (b) Comparison between time domain solution and steady solution.

4.3.6. $U=1.5{c_0}$

As shown in figure 15, for the moving body speed $U=1.5{c_0}$, the amplitude of the internal wave decreases further to $0.15h_{2}$. The time domain solution and the steady solution are also in good agreement in figure 15(b).

Figure 15. Generation of internal waves with the moving body speed $U=1.5c_{0}$, with $\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$. (a) Comparison of the internal wave profiles at different moments. (b) Comparison between time domain solution and steady solution.

Then we study the relationship between the moving body speed and the amplitude of the generated internal wave when the moving body speed increases from $U=1.242c_{0}$ to $U=1.5c_{0}$. As shown in figure 16, we observe that the generated internal wave amplitude decreases as the moving body speed increases. The relation, however, is nonlinear, unlike what was observed for $U\leq 1.24c_{0}$ (in figure 11).

Figure 16. Relationship between the generated internal wave amplitude $a_{1}$ and the moving body speed $U$, the moving speed changing from $U=1.242c_{0}$ to $U=1.5c_{0}$, with $\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$.

In figure 17, we combine figures 11 and 16. It is obvious that there is a critical speed between $U=1.241c_{0}$ and $U=1.242c_{0}$. The amplitude of the internal wave decreases significantly after exceeding the critical speed. For example, the amplitude of the internal wave under the moving body speed $U=1.241c_{0}$ reaches $1.23h_{2}$, but when the moving speed exceeds the speed $U=1.242c_{0}$, the amplitude of the generated internal wave decreases to $0.41h_{2}$ rapidly. This phenomenon was also observed by Grue et al. (1997) in the simulation of internal waves generated by a moving body on the free surface. It is interesting to note that similar observations were reported by Ertekin (1984), Ertekin, Webster & Wehausen (1986) and Ertekin, Qian & Wehausen (1990) for a surface ship and surface disturbance, and by Ertekin (1984) for a submerged bottom bump, although they were for a single layer of fluid.

Figure 17. Relationship between the generated internal wave amplitude $a_{1}$ and the moving body speed $U$, the moving speed changing from $U=0.8c_{0}$ to $U=1.5c_{0}$, with $\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$.

To assess how the dimensions of the bottom disturbance affect the amplitude of the internal wave generated and the critical speed, we consider three bottom disturbances with different major and minor axes. Since the amplitudes of the internal waves generated at different moving speeds are different, here we compare the amplitude of the leading wave generated when the moving speed is $U=1.0c_{0}$. By comparing the results of cases A and B in table 5, when the semi-major axis of the semi-ellipse is doubled, we find that the amplitude of the leading wave is slightly smaller, and the critical speed slightly increases. By comparing the results of cases A and C, when the semi-minor axis of the semi-ellipse is doubled, it is observed that the amplitude of the leading wave is significantly larger, and the critical speed increases.

Table 5. The effect of different dimensions of the bottom disturbance on the wave amplitude $a_{1}$ and the critical speed ${c_{{{critical}}}}$, where $a_{1}$ is the amplitude of the leading wave when the moving speed is $U=1.0c_{0}$. Here, $L_{1/2}$ is the semi-major axis of the semi-ellipse and $B_{1/2}$ is its semi-minor axis.

4.4. Internal waves generated by an unsteady moving bottom disturbance

The MCC-RL model can be applied to simulate the internal waves generated by an unsteady bottom disturbance. Here, we perform a numerical study on the internal waves generated by a moving body with a variable speed. The parameters are the same as those given in § 4.2 except for the moving speed of the bottom disturbance. The moving speed at different moments is shown in figure 18.

Figure 18. The unsteady moving speed of the body at different moments.

As shown in figure 18, at the first stage, the body is accelerated at a constant rate from 0 to $U=1.1c_{0}$ in the first 50 s. At the second stage, it moves at a constant speed $U=1.1c_{0}$ for the next 200 s. At the last stage, it is decelerated at a constant rate from $U=1.1c_{0}$ to 0 in the last 50 s. The internal wave profiles generated due to this unsteady motion, obtained by the MCC-RL model, are shown in figure 19 for different moments.

Figure 19. The internal wave profiles generated by an unsteady moving bottom, with $\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$, for (a) $t=50$ s, (b) $t=100$ s, (c) $t=200$ s, (d) $t=250$ s, (e) $t=300$ s and (f) $t=400$ s.

As shown in figure 19, at the end of the first stage ($t=50$ s), the internal solitary wave is not fully generated yet. At the second stage, several internal solitary waves are generated. At $t = 200$ s, it can be seen that the amplitudes of the third and fourth internal waves are relatively small, while at $t = 250$ s, the amplitude of the third wave increases. Meanwhile, the fourth and fifth internal solitary waves with smaller amplitudes are formed. At the last stage ($t = 300$ s), an obvious disturbance is observed right above the body, and the amplitude of the sixth wave is relatively small. After $t = 300$ s, the body is quiescent and the generated internal waves can propagate steadily as shown in $t = 400$ s.

5. Conclusions

In this paper, we apply the MCC-RL model to solve the time-varying bottom problems. The equations of the MCC-RL model with time-varying bottom are introduced, and the numerical algorithm for the time domain simulations is given. By testing the steady propagation of a large-amplitude internal solitary wave on the flat bottom, the accuracy of the time domain results and the numerical algorithm are verified.

We focus on the numerical simulations on the internal waves generated by a moving semi-ellipse (semi-major axis $10h_{2}$, semi-minor axis $0.1h_{2}$) on the bottom in a two-layer system ($\rho _1/\rho _2=0.7873$ and $h_1/h_2=4/1$); the conclusions are outlined below.

1. (i) The results of the MCC-RL model and Euler's solutions by Grue et al. (1997) match well in general for the case $U=1.1c_0$. Good agreements are also found between the time domain solution and the steady solution, which indicates that the generated internal waves are internal solitary waves indeed.

2. (ii) By changing the moving speed from $U=0.8c_{0}$ to $U=1.5c_{0}$, we find that there exists a critical speed between $U=1.241c_0$ and $U=1.242c_0$. There are significant differences on the generated internal waves when the speed is smaller than the critical speed and when the speed is greater than the critical speed.

3. (iii) When the moving speed is smaller than the critical speed, the internal solitary waves can be generated continuously, and they are in good agreement with the steady solutions of the internal solitary wave with flat bottom. With the moving body speed increasing from $U=0.8c_{0}$ to $U=1.241c_{0}$, the amplitudes of the generated internal solitary waves increase monotonically, and the wave speed is closer to the moving body speed.

4. (iv) When the moving speed is greater than the critical speed, only one internal wave right above the body is generated, whose speed is the same as the moving speed. Good agreement is found between the time domain solution and the steady solution of the internal wave with time-varying bottom. After exceeding the critical speed, the amplitude of the internal wave decreases significantly. When the moving body speed increases from $U=1.242c_{0}$ to $U=1.5c_{0}$, the amplitude of the generated internal waves decreases further, but nonlinearly.