1. Introduction
Internal solitary waves exist widely in various seas all over the world; see e.g. Jackson (Reference Jackson2007). Internal solitary waves are seen frequently in the northern part of the South China Sea, where Huang et al. (Reference Huang, Chen, Zhao, Zhang, Zhou, Yang and Tian2016) 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 (Reference Benjamin1966) and Miles (Reference Miles1980). 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. Reference Grue, Jensen, Rusas and Sveen1999; Kodaira et al. Reference Kodaira, Waseda, Miyata and Choi2016). Thus the KdV model is not suitable for describing large-amplitude internal waves due to the assumption of weak nonlinearity (Ostrovsky & Stepanyants Reference Ostrovsky and Stepanyants2005; Helfrich & Melville Reference Helfrich and Melville2006).
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 (Reference Miyata1985, Reference Miyata1988) and Choi & Camassa (Reference Choi and Camassa1999) (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 (Reference Choi and Camassa1996) 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 Reference Choi and Camassa1999; Camassa et al. Reference Camassa, Choi, Michallet, Rusas and Sveen2006; Barros & Gavrilyuk Reference Barros and Gavrilyuk2007; Jo & Choi Reference Jo and Choi2008; Kodaira et al. Reference Kodaira, Waseda, Miyata and Choi2016; la Forgia & Sciortino Reference la Forgia and Sciortino2019). Recently, Barros et al. (Reference Barros, Choi and Milewski2020) 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 (Reference la Forgia and Sciortino2020, Reference la Forgia and Sciortino2021) 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 (Reference Choi2022) 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 (Reference Choi2022) and the Euler's solutions.
For the space-varying bottom problems, Jo & Choi (Reference Jo and Choi2002) 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. (Reference Choi, Zhi and Barros2020) 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 (Reference Choi2000). 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 Reference Mercier, Vasseur and Dauxois2011; Duchene Reference Duchene2011) or by a seafloor disturbance, such as underwater landslides (see e.g. Brizuela, Filonov & Alford Reference Brizuela, Filonov and Alford2019). 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. (Reference Grue, Friis, Palm and Rusas1997) 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
$$\begin{gather} {u_{1,x}} + {w_{1,z}} = 0, \end{gather}$$
$$\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
$$\begin{gather} {u_{1,t}} + {u_1}{u_{1,x}} + {w_1}{u_{1,z}} = {{ - {p_{1,x}}} / {{\rho _1}}}, \end{gather}$$
$$\begin{gather}{u_{2,t}} + {u_2}{u_{2,x}} + {w_2}{u_{2,z}} = {{ - {p_{2,x}}} / {{\rho _2}}}, \end{gather}$$
$$\begin{gather}{w_{1,t}} + {u_1}{w_{1,x}} + {w_1}{w_{1,z}} = {{ - {p_{1,z}}} / {{\rho _1}}} - g, \end{gather}$$
$$\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
$$\begin{gather} {w_1} = 0,\quad{\rm at}\ z = {h_1}, \end{gather}$$
$$\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
$$\begin{gather} {w_2} = {\zeta_{,t}} + {u_2}{\zeta_{,x}},\quad {\rm at}\ z = {\zeta}(x,t), \end{gather}$$
$$\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
\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 (Reference Choi and Camassa1999) can be written in terms of four unknowns (
$\zeta,{\bar u_1},{\bar u_{2}},P$) as
$$\begin{gather} {\eta _{1,t}} + {\left( {{\eta _1}{{\bar u}_1}} \right)_{,x}} = 0, \quad{\eta _1} = {h_1} - \zeta, \end{gather}$$
$$\begin{gather}{\eta _{2,t}} + {\left( {{\eta _2}{{\bar u}_2}} \right)_{,x}} = 0, \quad {\eta _2} = {h_2} + \zeta, \end{gather}$$
$$\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}$$
$$\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
$$\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}$$
$$\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
\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 (Reference Choi2000), we obtain the equations for the MCC-RL model for the four unknowns (
$\zeta,{\bar u_1},{\bar u_2},P$) as
$$\begin{gather} {\eta _{1,t}} + {\left( {{\eta _1}{{\bar u}_1}} \right)_{,x}} = 0, \quad {\eta _1} = {h_1} - \zeta, \end{gather}$$
$$\begin{gather}{\eta _{2,t}} + {\left( {{\eta _2}{{\bar u}_2}} \right)_{,x}} = 0, \quad {\eta _2} = {h_2} + \zeta - b, \end{gather}$$
$$\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}$$
$$\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
\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
\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 (Reference Jo and Choi2008) 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
$$\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
\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
\begin{align} A &=-\frac{1}{3}\,{\eta _2}\left( {{\rho _1}{\eta _1} + {\rho _2}{\eta _2}} \right), \end{align}
\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}
\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}
\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:
$$\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}$$
$$\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
\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
$$\begin{gather} {\tilde A^{(i)}} =- {A^{(i)}}\,\frac{1}{{12{{\left( {\Delta x} \right)}^2}}} + {B^{(i)}}\,\frac{1}{{12\Delta x}}, \end{gather}$$
$$\begin{gather}{\tilde B^{(i)}} = { A ^{(i)}}\,\frac{{4}}{{3{{\left( {\Delta x} \right)}^2}}} - { B ^{(i)}}\,\frac{2}{{3\Delta x}}, \end{gather}$$
$$\begin{gather}{\tilde C^{(i)}} =- {A ^{(i)}}\,\frac{{5}}{{2{{\left( {\Delta x} \right)}^2}}} + {C ^{(i)}}, \end{gather}$$
$$\begin{gather}{\tilde D^{(i)}} = {A ^{(i)}}\,\frac{{4}}{{3{{\left( {\Delta x} \right)}^2}}} + { B ^{(i)}}\,\frac{2}{{3\Delta x}}, \end{gather}$$
$$\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 (Reference Zhao, Duan and Ertekin2014). 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 (Reference Jo and Choi2002). 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 (Reference Jo and Choi2008) and Choi et al. (Reference Choi, Zhi and Barros2020). On the other hand, the internal wave model can be regularized to eliminate shear instability; see e.g. Choi, Barros & Jo (Reference Choi, Barros and Jo2009), Lannes & Ming (Reference Lannes and Ming2015) and Duchêne, Israwi & Talhouk (Reference Duchêne, Israwi and Talhouk2016).
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
$$\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}$$
$$\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 Reference Longuet-Higgins and Cokelet1976; Fuhrman, Madsen & Bingham Reference Fuhrman, Madsen and Bingham2006; Zhao et al. Reference Zhao, Duan and Ertekin2014).
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 (Reference Jo and Choi2002). 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. (Reference Grue, Jensen, Rusas and Sveen1999). 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. (Reference Grue, Friis, Palm and Rusas1997). 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. (Reference Grue, Jensen, Rusas and Sveen1999), 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 (Reference Choi and Camassa1999).
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. (Reference Grue, Jensen, Rusas and Sveen1999) 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. Reference Camassa, Choi, Michallet, Rusas and Sveen2006). In figure 4(b), 
$c_0$ is the linear long-wave speed, which is defined as
\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 Reference Wang2019). 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 Reference Grimshaw2010; Grimshaw & Maleewong Reference Grimshaw and Maleewong2016) and in a continuously density-stratified fluid (Grimshaw & Smyth Reference Grimshaw and Smyth1986; Grimshaw, Chan & Chow Reference Grimshaw, Chan and Chow2002; Grimshaw & Helfrich Reference Grimshaw and Helfrich2018). 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 Reference Helfrich and Melville2006).
In the simulation, we use the same parameters given by Grue et al. (Reference Grue, Friis, Palm and Rusas1997), 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. (Reference Grue, Jensen, Rusas and Sveen1997) 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. (Reference Grue, Friis, Palm and Rusas1997), 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 (Reference Choi and Camassa1999), as are also shown below:
$$\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}$$
$$\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. (Reference Grue, Friis, Palm and Rusas1997).
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 (Reference Melville and Helfrich1987), 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
$$\begin{gather} \zeta (x,t) = \zeta (X), \end{gather}$$
$$\begin{gather}b(x,t) = b(X), \end{gather}$$
$$\begin{gather}{\bar u_1}(x,t) = {\bar u_1}(X), \end{gather}$$
$$\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
\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. (Reference Zhao, Ertekin, Duan and Webster2016) and Duan et al. (Reference Duan, Wang, Zhao, Ertekin and Kim2018). 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. (Reference Grue, Friis, Palm and Rusas1997) 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 (Reference Ertekin1984), Ertekin, Webster & Wehausen (Reference Ertekin, Webster and Wehausen1986) and Ertekin, Qian & Wehausen (Reference Ertekin, Qian and Wehausen1990) for a surface ship and surface disturbance, and by Ertekin (Reference Ertekin1984) 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.
(i) The results of the MCC-RL model and Euler's solutions by Grue et al. (Reference Grue, Friis, Palm and Rusas1997) 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.(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.(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.(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.
Acknowledgements
The authors are grateful to the anonymous referees for their comments that improved our paper.
Funding
The work of B.Z., W.D. and Z.W. is supported by the National Natural Science Foundation of China (no. 12202114), the China Postdoctoral Science Foundation (no. 2022M710932), the State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology (no. LP2202), the Fundamental Research Funds for the Central Universities (no. 3072022FSC0101), the Qingdao Postdoctoral Application Project and the Heilongjiang Touyan Innovation Team Programme.
Declaration of interests
The authors report no conflict of interest.
