1. Introduction
Wave dynamics in the oceanic environment with density stratified structures, known as internal waves, are essential for untangling many physical and environmental processes in oceanography. The topographic effects on generation, propagation and interactions of internal waves of different modes have attracted great interest in the past few decades. The corresponding mechanisms are vital in physical oceanography and coastal engineering applications. The bottom boundary, with varying forms in the typical oceanic area such as strait, submerged reef and seacoast, represents a specific and primary source of disturbance on internal waves. Though it will increase complexity, the necessity of involving real bathymetry in the internal wave system has inspired current research.
Numerically solving the classical water-wave problem, i.e. the full Euler equations composed of the Laplace equation governing the motion of the bulk fluid, the impermeability condition at the bottom and the kinematic and dynamic conditions on the free surface, is usually challenging. Its core difficulty lies in accurately obtaining the normal velocity on the free surface at every time step. In a two-dimensional fluid domain (corresponding to a one-dimensional free surface), numerical methods designed based on complex analysis in one variable perform well, with or without the inclusion of bottom topography. However, other types of numerical methods that accommodate more general cases have to be developed since the complex function theory cannot be applied or generalized to three-dimensional problems. A large category of methods is to reformulate the system with the surface variables using a Dirichlet–Neumann map for the velocity potential so as to reduce the dimension of the problem. The advantage of this type of method is avoiding the expensive numerical solution of the Laplace equation in an irregular time-varying domain, and the key to such methods lies in the accurate and efficient computation of the Dirichlet-to-Neumann (DtN) operator that maps the given Dirichlet data to the associate Neumann data at the free surface. Wilkening & Vasan (Reference Wilkening and Vasan2015) summarized four existing forms of the DtN operator and compared them with the boundary integral collocation method from some numerical perspectives.
The first characterization of the DtN operator is the Craig–Sulem expansion pioneered by Craig & Selum (Reference Craig and Selum1993). They proposed expressing the DtN operator as a Taylor series by assuming the analyticity of the wave profile and using a pseudo-spectral method for numerical computations. Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005) and Guyenne & Nicholls (Reference Guyenne and Nicholls2005) extended the Craig–Sulem expansion to water waves over a periodically varying bottom, showing that the topographic effect introduces an additional term in the DtN operator. This additional operator also has a convergence Taylor series when the deformation in the bottom boundary is analytic and small. With the aid of this extended DtN operator, Craig et al. (Reference Craig, Guyenne, Nicholls and Sulem2005) presented, in various asymptotic limits, a systematic derivation of long-wave nonlinear models with a non-trivial bathymetry based on the Hamiltonian perturbation analysis. Guyenne & Nicholls (Reference Guyenne and Nicholls2005) implemented the explicit recursion formula of the topographic operator to simulate the evolution of solitary waves over plane slopes. They noticed the instability at high wavenumbers, presumably owing to aliasing errors or the ill-conditioning of the operator expansion. The Craig–Sulem expansion of the DtN operator was also generalized to a moving bottom topography by Guyenne & Nicholls (Reference Guyenne and Nicholls2008) to mimic the earthquake- and landslide-generated tsunamis. The authors rewrote the formal expression of Taylor series terms to eliminate the explicit appearance of extremely unstable operators. Upon noticing that the topographic operator in the works mentioned above (Craig et al. Reference Craig, Guyenne, Nicholls and Sulem2005; Guyenne & Nicholls Reference Guyenne and Nicholls2005, Reference Guyenne and Nicholls2008) already includes a smoothing operator, Cathala (Reference Cathala2016) derived a non-local shallow-water model in which only smoothing contribution is involved for flows over a non-smooth topography and showed numerically a satisfactory comparison with that of Nachbin (Reference Nachbin2003) for the case of a polygonal topography.
The second and third methods, the Ablowitz–Fokas–Musslimani (AFM) implicit formulation and the dual AFM (AFM*) formulation, for characterizing the DtN operator, are homologous. To construct a global relation between the Dirichlet and Neumann data, basic hyperbolic functions and Green's second identity are used to obtain integral equations. Solving the implicit or explicit non-local characterization is equivalent to summing all terms in the Craig–Sulem expansion. Milewski (Reference Milewski1998) proposed an early derivation with weakly nonlinear approximation and asymptotic expansion of velocity potential, where a Fourier integral equation was obtained as an alternative formulation of the water-wave problem involving topographic and capillary effects. From the mathematical point of view, Fokas (Reference Fokas2000) introduced a general approach for boundary value problems for linear and integrable nonlinear partial differential equations and showed that the last step of this approach was to analyse the global relation between the solution and its derivatives on the boundary, which could be obtained by constructing and integrating Green's function. Based on previous works, Ablowitz, Fokas & Musslimani (Reference Ablowitz, Fokas and Musslimani2006) presented the AFM non-local formulation of water waves with non-local spectral equations where both surface tension and bottom topography were taken into consideration. The AFM* formulation is a dual version of the AFM method (noticing the adjoint property of the DtN operator) initially derived by Ablowitz & Haut (Reference Ablowitz and Haut2008) for interfacial waves in the classic two-layer system. After that, considerable effort has been devoted to utilizing/generalizing the AFM* method from multiple aspects. Of note is the work of Haut & Ablowitz (Reference Haut and Ablowitz2009) who derived the Benny–Luke and intermediate long-wave equations for interfacial waves with a free surface and numerically investigated lump-type solutions, Vasan & Deconinck (Reference Vasan and Deconinck2013) who used the fully nonlinear AFM* formulation to study the inverse problem for bathymetry detection in a two-dimensional system, and Andrade & Nachbin (Reference Andrade and Nachbin2018) who developed a three-dimensional DtN operator based upon the linearized AFM* formulation. It is worth mentioning that in the linear approximation, the Craig-Sulem expansion and AFM formulation give the same result for the surface water-wave problem, regardless of the amplitude of bottom topography (see Craig et al. Reference Craig, Guyenne, Nicholls and Sulem2005).
The last characterization of the DtN operator is the transformed field expansion method. It is, in fact, an application of boundary-flattening change of variables to avoid the amplification of a round-off error in computations of the Craig–Sulem expansion terms. Nicholls & Reitich (Reference Nicholls and Reitich2001, Reference Nicholls and Reitich2006) considered a three-dimensional problem and proved the change of variables can lead to an accurate and well-conditioned perturbative calculation. These works do not include the topographic effect.
Except for the DtN operator, there also exist other reformulations of the water-wave problem in the spirit of Zakharov's equations, among which the most famous one is the high-order spectral (HOS) method pioneered by Dommermuth & Yue (Reference Dommermuth and Yue1986). The basis of the HOS method is the mode coupling idea, and the topographic effect can also be incorporated into the coupling. Instead of introducing the DtN operator, a series of boundary value problems are given at each order of the Taylor expansion of the velocity potential 
$\varPhi$. At 
$m$th order, with the mode coupling approach, 
$\varPhi ^{(m)}$ has an eigenfunction expansion of free modes. The time advancement of the amplitude of each mode is computed with free-surface conditions. The advantage of this method is the accessibility to the evolution of each mode calculated. However, the restriction on the system is unavoidable; for example, the bottom deformation and free-surface displacement are demanded to be small. Liu & Yue (Reference Liu and Yue1998) performed the numerical simulation with the HOS method to investigate the generalized Bragg scattering in three-dimensional cases, and Alam, Liu & Yue (Reference Alam, Liu and Yue2009a,Reference Alam, Liu and Yueb) extended this method for a two-fluid system with a free surface for analyses of the Bragg resonance over bottom ripples.
When it comes to a more realistic oceanic situation, the complexity of the internal wave problem over real bathymetry (which could be smoother, slowly varying or periodic) poses more challenges. To the best of the authors’ knowledge, among the above methods, only AFM* and HOS have been generalized to the direct numerical simulation of internal waves. Most existing studies focus on deriving the reduced models, such as the Korteweg–de Vries (KdV) and Kadomtsev–Petviashvili (KP) equations. For models considering seabed terrain, Chen & Liu (Reference Chen and Liu1996) derived the modified KdV equation for a two-fluid system with a critical depth ratio and investigated the interfacial waves over random topography. They found that the randomness of the topography with white Gaussian noise distribution results in an averaged solitary wave, which gradually deforms into a spreading Gaussian wavepacket. Lynett & Liu (Reference Lynett and Liu2002) considered a two-layer system with a free surface and proposed a weakly dispersive, weakly nonlinear and depth-integrated model for interfacial waves passing over real bathymetry. However, even though the authors did not use the rigid-lid approximation, the amplitude of waves on the free surface was much smaller than internal waves due to a small density difference between two layers; thus, only internal waves were numerically solved and compared with satellite images for model verification. Nonlinear models for internal waves in a continuously stratified density environment were introduced by Grimshaw (Reference Grimshaw1981). Based on generalized Lagrangian variables, Grimshaw separated the motion of the fluid into a background shear flow and a relative perturbation. With asymptotic analysis and method of multiple scales, the variable coefficient KdV and Benjamin–Davis–Ono equations were derived for shallow and deep water cases, respectively. Later on, Helfrich & Melville (Reference Helfrich and Melville1986) extended the model to include bottom topography and its application to tidal flows can be referred to Griffiths & Grimshaw (Reference Griffiths and Grimshaw2007). Yuan et al. (Reference Yuan, Grimshaw, Johnson and Chen2018b) derived a variable coefficient KP equation in a three-dimensional system and compared the evolution of internal solitary waves with those of the Massachusetts Institute of Technology general circulation model. Grimshaw's method can also be generalized to study multimode internal waves. The variable coefficient KdV equation was used to simulate the evolution of a mode 2 internal solitary wave over a slope-shelf topography by Yuan, Grimshaw & Johnson (Reference Yuan, Grimshaw and Johnson2018a). Moreover, Liu, Grimshaw & Johnson (Reference Liu, Grimshaw and Johnson2019) used a limited linear coupling model between mode 1 and mode 2 to analyse the mode-2 generation as mode 1 propagates shoreward.
The oceanic environment usually shows pycnocline, the density of which varies rapidly; thus, the ocean can approximately be separated into three layers with different typical densities. The stratification is continuous in the pycnocline, with infinitely many interfaces presented so that multimode internal waves exist. Consequently, the three-layer fluid system with two interfaces may be the simplest idealization to be analogous to the real situation and is considered in this paper with the further assumption of being incompressible, inviscid and irrotational for each layer fluid; thus, interactions between mode 1, mode 2 and bathymetry can be analysed. When the density of the top layer is negligible (for example, the air–water interface), the three-layer system reduces to the two-fluid system with a free surface on the top, which can be considered a limiting case of the three-layer configuration. Multimode nonlinear internal waves in the three-layer configuration with flat boundaries have attracted considerable interest in recent years, and particular attention has been paid to the mode 2 solitary waves without oscillatory tails. For example, mode 2 free solitary waves have been found in the strongly nonlinear Miyata–Choi–Camassa type equations by Barros, Choi & Milewski (Reference Barros, Choi and Milewski2020) for the shallow-water scenario and by Wang, Wang & Yuan (Reference Wang, Wang and Yuan2022) for fluids of great depth, and solutions in the former situation have also been confirmed in the full Euler equations by Doak, Barros & Milewski (Reference Doak, Barros and Milewski2022). In the present work, the DtN operator method is generalized to a three-layer system. The ansatz for the topographic component follows the expression of Milewski (Reference Milewski1998). Inspired by Andrade & Nachbin (Reference Andrade and Nachbin2018), internal waves are assumed to be linear, but no restriction is imposed on the topography. Both two- and three-dimensional evolutions of multimode internal waves over bathymetry are simulated, and a Galerkin approximation is applied to reduce the computational cost when calculating topographic coefficients.
The rest of the paper is structured as follows. In § 2 the mathematical formulation of interfacial waves in a three-layer fluid system passing over bottom topography, together with the linear dispersion relation, is presented. The numerical method based on a generalization of Andrade & Nachbin (Reference Andrade and Nachbin2018) by introducing two DtN operators is constructed in § 3. Corresponding lateral boundary conditions and initialization are also explained. Section 4 presents the main simulation results for multimode internal waves in both two and three dimensions. Finally, the concluding remarks are given in § 5.
2. Mathematical formulation
A system composed of three superposed inviscid, immiscible and incompressible fluids separated by two interfaces is considered in three dimensions (see figure 1 for a sketch). The system is bounded above by a flat rigid wall and below by variable topography. The fluid in each layer, which is assumed to be irrotational, has finite depth and constant density. A Cartesian coordinate system 
$(x,y,z)$ is introduced such that 
$x$ and 
$y$ are horizontal coordinates and the 
$z$ axis is in the opposite direction of gravity. Properties in the lower layer are designated with subscript 1 and those in the middle and upper layers are analogously assigned with subscripts 2 and 3, respectively. The depth of the lower layer, 
$H_1$, is defined with the vertical distance between the position of the lower interface at rest and the reference bottom position 
$z=0$. The bottom topography is expressed by 
$z=b(x,y)$ and the two interfaces are designated by 
$z=H_1+\eta _1(x,y,t)$ and 
$z=H_1+H_2+\eta _2(x,y,t)$, with 
$t$ being the time. The lighter fluid is always placed on top of the heavier one to avoid the Rayleigh–Taylor instability, namely 
$\rho _3<\rho _2<\rho _1$ in figure 1.
Figure 1. Sketch of interfacial waves propagating over large bottom topography in a three-layer density stratified system.
Applying the potential theory, the motions of fluid bodies are governed by the Laplace equations, which can be written as
\begin{equation} \left.\begin{array}{c} \phi_{1xx}+\phi_{1yy}+\phi_{1zz}=0 \quad \text{for } b(x,y)< z< H_1+\eta_1(x,y,t),\\ \phi_{2xx}+\phi_{2yy}+\phi_{2zz}=0 \quad \text{for } H_1+\eta_1(x,y,t)< z< H_1+H_2+\eta_2(x,y,t),\\ \phi_{3xx}+\phi_{3yy}+\phi_{3zz}=0 \quad \text{for } H_1+H_2+\eta_2(x,y,t)< z< H_1+H_2+H_3. \end{array}\right\} \end{equation}The impermeability boundary conditions on the top and bottom walls can be expressed as
\begin{equation} \left.\begin{array}{c} \phi_{1z}-\phi_{1x}b_x-\phi_{1y}b_y=0 \quad \text{on } z=b(x,y),\\ \phi_{3z}=0 \quad \text{on } z=H_1+H_2+H_3. \end{array}\right\} \end{equation}
On the interfaces, the kinematic boundary conditions governing the time evolutions of the wave displacements, 
$\eta _1(x,y,t)$ and 
$\eta _2(x,y,t)$, read
\begin{equation} \left.\begin{array}{c} \eta_{1t}+\eta_{1x}\phi_{1x}+\eta_{1y}\phi_{1y}-\phi_{1z}=0 \quad \text{on } z=H_1+\eta_1(x,y,t), \\ \eta_{1t}+\eta_{1x}\phi_{2x}+\eta_{1y}\phi_{2y}-\phi_{2z}=0 \quad \text{on } z=H_1+\eta_1(x,y,t), \\ \eta_{2t}+\eta_{2x}\phi_{2x}+\eta_{2y}\phi_{2y}-\phi_{2z}=0 \quad \text{on } z=H_1+H_2+\eta_2(x,y,t), \\ \eta_{2t}+\eta_{2x}\phi_{3x}+\eta_{2y}\phi_{3y}-\phi_{3z}=0 \quad \text{on } z=H_1+H_2+\eta_2(x,y,t). \end{array}\right\} \end{equation}Finally, the Bernoulli principle yields the dynamic conditions on the interfaces, namely,
\begin{equation} \left.\begin{array}{c} \begin{aligned} & \phi_{1t}-R_1\phi_{2t}+g(1-R_1)\eta_1+ \tfrac{1}{2}(\phi_{1x}^2+\phi_{1y}^2+\phi_{1z}^2-R_1\phi_{2x}^2 -R_1\phi_{2y}^2-R_1\phi_{2z}^2)=0 \\ & \quad\text{on } z=H_1+\eta_1(x,y,t), \\ & \phi_{2t}-R_2\phi_{3t}+g(1-R_2)\eta_2+\tfrac{1}{2} (\phi_{2x}^2+\phi_{2y}^2+\phi_{2z}^2-R_2\phi_{3x}^2-R_2\phi_{3y}^2-R_2\phi_{3z}^2) =0 \\ & \quad\text{on } z=H_1+H_2+\eta_2(x,y,t), \end{aligned} \end{array}\right\} \end{equation}
where 
$R_1=\rho _2/\rho _1<1$ and 
$R_2=\rho _3/\rho _2<1$. We remark that though there is no smallness assumption on topographic relief, a restriction, 
$b(x,y)< H_1+\eta _1$, is still imposed in the present problem.
2.1. Non-dimensional form and linear model
We non-dimensionalize the system using the scalings
\begin{align} x'=lx,\quad y'=ly,\quad z'=H_1z,\quad \eta'=a\eta,\quad b'=H_1b,\quad t'=\frac{l}{\sqrt{gH_1}}t,\quad \phi'=\frac{agl}{\sqrt{gH_1}}\phi, \end{align}
where 
$l$ is the reference wavelength and 
$a$ is the typical wave amplitude. The small dimensionless parameter 
$\epsilon =a/H_1\ll 1$ arises along with 
$\mu =H_1/l$ for which, notably, no assumption is imposed (which is different from the Boussinesq scaling). Here 
$\mu$ will be used as a parameter whose value is preassigned, and we choose 
$\mu =0.1$ in most computations in the present paper.
Dropping the apostrophes of dimensionless variables, the governing equations and impermeability boundary conditions in non-dimensional form can be recast to
\begin{equation} \left.\begin{array}{c} \mu^2(\phi_{1xx}+\phi_{1yy})+\phi_{1zz}=0 \quad \text{for } b< z<1+\epsilon\eta_1,\\ \mu^2(\phi_{2xx}+\phi_{2yy})+\phi_{2zz}=0 \quad \text{for } 1+\epsilon\eta_1< z<1+\varTheta_1+\epsilon\eta_2,\\ \mu^2(\phi_{3xx}+\phi_{3yy})+\phi_{3zz}=0 \quad \text{for } 1+\varTheta_1+\epsilon\eta_2< z<1+\varTheta_1+\varTheta_2,\\ \phi_{1z}-\mu^2(\phi_{1x}b_x+\phi_{1y}b_y)=0 \quad \text{on } z=b,\\ \phi_{3z}=0 \quad \text{on } z=1+\varTheta_1+\varTheta_2, \end{array}\right\} \end{equation}
where 
$\varTheta _1=H_2/H_1$ and 
$\varTheta _2=H_3/H_1$ are the depth ratios. The kinematic and dynamic boundary conditions on the interfaces are, on 
$z=1+\epsilon \eta _1$,
\begin{equation} \left.\begin{aligned} &\displaystyle \eta_{1t}=\dfrac{1}{\mu^2}\phi_{1z}-\epsilon(\phi_{1x}\eta_{1x}+\phi_{1y}\eta_{1y})=\dfrac{1}{\mu^2}\phi_{2z}-\epsilon(\phi_{2x}\eta_{1x}+\phi_{2y}\eta_{1y}),\\ &\displaystyle \phi_{1t}-R_1\phi_{2t}+(1-R_1)\eta_1\\ &\quad +\dfrac{\epsilon}{2} \left(\phi_{1x}^2+\phi_{1y}^2+\dfrac{1}{\mu^2}\phi_{1z}^2-R_1\phi_{2x}^2-R_1\phi_{2y}^2-\dfrac{R_1}{\mu^2}\phi_{2z}^2\right)=0, \end{aligned}\right\} \end{equation}
and on 
$z=1+\varTheta _1+\epsilon \eta _2$,
\begin{equation} \left.\begin{aligned} &\displaystyle \eta_{2t}=\dfrac{1}{\mu^2}\phi_{2z}-\epsilon(\phi_{2x}\eta_{2x}+\phi_{2y}\eta_{2y})=\dfrac{1}{\mu^2}\phi_{3z}-\epsilon(\phi_{3x}\eta_{2x}+\phi_{3y}\eta_{2y}),\\ &\displaystyle \phi_{2t}-R_2\phi_{3t}+(1-R_2)\eta_2\\ &\quad +\dfrac{\epsilon}{2}\left(\phi_{2x}^2+\phi_{2y}^2+\dfrac{1}{\mu^2}\phi_{2z}^2-R_2\phi_{3x}^2-R_2\phi_{3y}^2-\dfrac{R_2}{\mu^2}\phi_{3z}^2\right)=0. \end{aligned}\right\} \end{equation}
As the linear dispersive system has been proven to be well suited for the investigation of small-amplitude waves in a realistic physical scale, following the previous work (Andrade & Nachbin Reference Andrade and Nachbin2018), the quadratic terms of 
$O(\epsilon )$ in (2.6)–(2.8) will be dropped. It is worth noting that the topographic effect can be retained in the bottom boundary condition by regarding 
$\mu$ as a variable parameter. Finally, the linearized system to be solved can be rewritten as
\begin{gather} \left.\begin{array}{c} \mu^2(\phi_{1xx}+\phi_{1yy})+\phi_{1zz}=0 \quad \text{for } b< z<1,\\ \mu^2(\phi_{2xx}+\phi_{2yy})+\phi_{2zz}=0 \quad \text{for } 1< z<1+\varTheta_1,\\ \mu^2(\phi_{3xx}+\phi_{3yy})+\phi_{3zz}=0 \quad \text{for } 1+\varTheta_1< z<1+\varTheta_1+\varTheta_2, \end{array}\right\} \end{gather}
\begin{gather} \left.\begin{array}{c} \phi_{1z}-\mu^2(\phi_{1x}b_x+\phi_{1y}b_y)=0 \quad \text{on } z=b,\\ \phi_{3z}=0 \quad \text{on } z=1+\varTheta_1+\varTheta_2, \end{array}\right\} \end{gather}
\begin{gather} \left.\begin{array}{c} \displaystyle \eta_{1t}-\dfrac{1}{\mu^2}\phi_{1z}=\eta_{1t}-\dfrac{1}{\mu^2}\phi_{2z}=0 \quad \text{on } z=1,\\ \displaystyle \eta_{2t}-\dfrac{1}{\mu^2}\phi_{2z}=\eta_{2t}-\dfrac{1}{\mu^2}\phi_{3z}=0 \quad \text{on } z=1+\varTheta_1, \end{array}\right\} \end{gather}
\begin{gather} \left.\begin{array}{c} \phi_{1t}-R_1\phi_{2t}+(1-R_1)\eta_1=0 \quad \text{on } z=1,\\ \phi_{2t}-R_2\phi_{3t}+(1-R_2)\eta_2=0 \quad \text{on } z=1+\varTheta_1. \end{array}\right\} \end{gather}2.2. Dispersion relation
We investigate the characteristics of the linear system with trivial bottom topography, i.e. 
$b(x,y)=0$. Since the system is isotropic in horizontal directions, it suffices to focus on the two-dimensional problem. Considering monochromatic incident waves, the wavenumber is related to the frequency through the dispersion relation. Assuming that two interfacial waves with sinusoidal profiles propagate in the 
$x$ direction from left to right without phase difference, namely
\begin{equation} \eta_1=a_1\sin{(kx-\omega_0t)},\quad \eta_2=a_2\sin{(kx-\omega_0t)}, \end{equation}
we then calculate the dispersion relation and the expression of 
$a_2/a_1$ as a function of wavenumber and angular velocity 
$(k,\omega _0)$. First, for a flat bottom, the impermeability boundary conditions in (2.10) become
\begin{equation} \left.\begin{array}{c} \phi_{1z}=0 \quad \text{on } z=0, \\ \phi_{3z}=0 \quad \text{on } z=1+\varTheta_1+\varTheta_2. \end{array}\right\} \end{equation}Substituting (2.13a,b) into (2.11), the vertical derivatives of velocity potentials on two interfaces can be obtained as
\begin{equation} \left.\begin{array}{c} \phi_{1z}=\phi_{2z}={-}\mu^2\omega_0a_1\cos{(kx-\omega_0t)} \quad \text{on } z=1,\\ \phi_{2z}=\phi_{3z}={-}\mu^2\omega_0a_2\cos{(kx-\omega_0t)} \quad \text{on } z=1+\varTheta_1. \end{array}\right\} \end{equation}Then, the Laplace equations in (2.9) can be solved, and the velocity potentials can be expressed as
\begin{equation} \left.\begin{array}{c}
\displaystyle \phi_1={-}\dfrac{\mu\omega_0a_1}{k\sinh{(\mu
k)}}\cos{(kx-\omega_0t)}\cosh{(\mu kz)},\\ \begin{aligned}
\phi_2 & =\dfrac{\mu\omega_0}{k\sinh^2{(\mu
k\varTheta_1)}}\cos{(kx-\omega_0t)}[(a_1-a_2\cosh{(\mu
k\varTheta_1)})\sinh{(\mu k(z-1))}\\ & \quad
+(a_2-a_1\cosh{(\mu k\varTheta_1)})\sinh{(\mu
k(z-1-\varTheta_1))}], \end{aligned}\\ \displaystyle
\phi_3=\dfrac{\mu\omega_0a_2}{k\sinh(\mu
k\varTheta_2)}\cos{(kx-\omega_0t)}\cosh{(\mu
k(z-1-\varTheta_1-\varTheta_2))}. \end{array}\right\}
\end{equation}Substituting the velocity potentials and wave profiles into the dynamic boundary conditions on two interfaces (see (2.12)), and after some algebra, one obtains
\begin{equation} \frac{a_2}{a_1}=\frac{\sinh(\mu k\varTheta_1)}{R_1\tanh(\mu k)} +\cosh(\mu k\varTheta_1)+\left(1-\frac{1}{R_1}\right) \frac{k}{\mu\omega_0^2}\sinh(\mu k\varTheta_1) \end{equation}and
\begin{equation} A\left(\frac{\mu\omega_0^2}{k^2}\right)^2+B\frac{\mu\omega_0^2}{k^2}+C=0, \end{equation}with
\begin{align} \left.\begin{array}{c} \displaystyle \dfrac{A}{k^2}=\dfrac{\tanh(\mu k\varTheta_2)+R_2\tanh(\mu k\varTheta_1)}{R_1}+\tanh(\mu k)[\tanh(\mu k\varTheta_1)\tanh(\mu k\varTheta_2)+R_2],\\ \begin{aligned} \dfrac{B}{k} & =\left(\dfrac{R_2}{R_1}-\dfrac{1}{R_1}\right) \tanh(\mu k\varTheta_1)\tanh(\mu k\varTheta_2) +\left(R_2-\dfrac{1}{R_1}\right)\tanh(\mu k)\tanh(\mu k\varTheta_2)\\ & \quad +\left(R_2-\dfrac{R_2}{R_1}\right)\tanh(\mu k)\tanh(\mu k\varTheta_1), \end{aligned}\\ \displaystyle C=(R_2-1)\left(1-\dfrac{1}{R_1}\right)\tanh(\mu k)\tanh(\mu k\varTheta_1)\tanh(\mu k\varTheta_2). \end{array}\right\} \end{align}
Equation (2.17) describes the ratio between wave displacements on two interfaces. Mode 1 internal waves will be obtained when 
$a_2/a_1$ is positive, and waves of mode 2 will be identified the other way around (i.e. 
$a_2/a_1<0$). The dispersion relation in a quadratic form is given in (2.18), and the unknown to be solved is 
$\mu \omega _0^2/k^2=\mu c_\phi ^2$, where 
$c_\phi$ is termed the phase velocity. With the fixed parameters 
$(R_1, R_2, \varTheta _1, \varTheta _2, \mu )$ and specific wavenumber 
$k$, two solutions can be obtained representing phase velocities for mode 1 and mode 2 internal waves, respectively. Some examples are given in figure 2.
Figure 2. Characteristics of monochromatic waves in three-layer configurations, including the dispersion relation between 
$k$ and 
$c_{\phi }$ and the ratio between wave displacements on two interfaces as a function of wavenumber: mode 1 internal waves (single lines) and mode 2 internal waves (lines with marks). Results are shown for (a,b) 
$(R_1, R_2, \mu )=(0.6,0.6,0.1)$ with 
$(\varTheta _1, \varTheta _2)=(1,1)$ (black solid lines), 
$(\varTheta _1, \varTheta _2)=(3,1)$ (blue dashed lines) and 
$(\varTheta _1, \varTheta _2)=(1,2)$ (red dotted lines); (c,d) 
$(R_1, \varTheta _1, \varTheta _2, \mu )=(0.6,1,1,0.1)$ with 
$R_2=0.6$ (black solid lines), 
$R_2=0.4$ (blue dashed lines) and 
$R_2=0.2$ (red dotted lines).
In figure 2(a) the curves of 
$c_{\phi }(k)$ for mode 1 and mode 2 approach each other when the depth of the middle layer increases. On the contrary, the distance between the two curves slightly increases when the depth of the upper layer increases. The density of the upper layer impacts the phase velocity for both modes in figure 2(c) (decreasing in 
$R_2$ leads to an increase in velocity for a fixed wavenumber), and the influence for mode 1 is more prominent. The ratio of wave displacement is monotone and can be ascending or descending. The opposite tendency of this ratio is observed in figures 2(b) and 2(d) for mode 1 and mode 2 waves. For a thick middle layer, mode 1 waves show increasing displacement ratio curves, while the ratio function decreases with wavenumber for a thick upper layer. The variable tendency can also be influenced by density ratios.
The solvability of the linear dispersion relation is worth attention. The solution, 
$\mu c_\phi ^2$, solved from (2.18) with pre-given system parameters should be non-negative; otherwise, the system becomes linearly ill-posed. It can be easily inferred from (2.19) that, for all 
$k>0$, 
$A>0$, 
$B<0$ and 
$C>0$. Assuming 
$B^2-4AC\geqslant 0$, the quadratic equation has real solutions giving rise to positive phase velocities. However, in the domains of the system parameters, i.e. 
$\varTheta _1,\varTheta _2\in (0,+\infty )$ and 
$R_1,R_2\in (0,1)$, it is not easy to prove that the discriminant is always non-negative. We consider the asymptotic limiting situations to elaborate on the potential restriction on system parameters. When 
$k\to +\infty$, it is not difficult to prove the positivity of the discriminant. We then consider the behaviour for a small wavenumber; asymptotic expansion with 
$k\to 0^+$ is performed for the coefficients in (2.19), which gives
\begin{equation} \left.\begin{array}{c} \displaystyle A\sim \mu k^3\left(\dfrac{\varTheta_2}{R_1}+\dfrac{\varTheta_1R_2}{R_1}+R_2\right),\\ \displaystyle B\sim \mu^2k^3\left[\varTheta_1\varTheta_2\left(\dfrac{R_2}{R_1}-\dfrac{1}{R_1}\right)+ \varTheta_2\left(R_2-\dfrac{1}{R_1}\right)+\varTheta_1R_2 \left(1-\dfrac{1}{R_1}\right)\right],\\ \displaystyle C\sim \mu^3k^3\varTheta_1\varTheta_2\left(R_2-1+\dfrac{1}{R_1}-\dfrac{R_2}{R_1}\right). \end{array} \right\} \end{equation}
By taking the limit 
$k\to 0^+$, the two solutions to the quadratic equation read
\begin{equation} \mu c_\phi^2=\mu\frac{\pm\sqrt{\varDelta}- \left[\varTheta_1\varTheta_2\left(\dfrac{R_2}{R_1}-\dfrac{1}{R_1}\right) +\varTheta_2\left(R_2-\dfrac{1}{R_1}\right)+\varTheta_1R_2 \left(1-\dfrac{1}{R_1}\right)\right]}{2\left(R_2+ \dfrac{\varTheta_1R_2}{R_1}+\dfrac{\varTheta_2}{R_1}\right)}, \end{equation}where
\begin{align} \varDelta&=\left[\varTheta_1\varTheta_2\left(\frac{R_2}{R_1}-\frac{1}{R_1}\right)+ \varTheta_2\left(R_2-\frac{1}{R_1}\right)+\varTheta_1R_2 \left(1-\frac{1}{R_1}\right)\right]^2 \nonumber\\ &\quad -4\varTheta_1\varTheta_2\left(R_2-1+\frac{1}{R_1}- \frac{R_2}{R_1}\right)\left(\frac{\varTheta_2}{R_1}+ \frac{\varTheta_1R_2}{R_1}+R_2\right). \end{align}
We can numerically show that, for small wavenumbers, the system remains feasible under most circumstances, and both modes of internal waves exist. Figure 3 illustrates the distribution of mode-2 phase velocities (which is the smaller solution) solved from (2.21) with different system parameters. Two particular parameters, for example, depth ratios in figure 3(a), are fixed, and the phase velocity is investigated in the plane of the other two parameters. Apparently, in figure 3 the phase velocity as 
$k\to 0^+$ is physical in all the areas considered. When the depth ratio tends to 0, or the density ratio tends to 1, 
$c_\phi ^2$ tends to 0. There is no indication that the three-layer system will become linearly unstable for a wide range of parameters for these two asymptotic limiting cases.
Figure 3. (a) Distribution of 
$c_\phi ^2$ in the 
$R_1$–
$R_2$ plane for 
$\varTheta _1=\varTheta _2=1$ and 
$\mu =0.1$. (b) Distribution of 
$c_\phi ^2$ in the 
$\varTheta _1$–
$\varTheta _2$ plane for 
$R_1=R_2=0.6$ and 
$\mu =0.1$. The feasibility region covers all the presented areas.
We finally remark that by taking the limit 
$\varTheta _2 \to +\infty$ and 
$R_2 \to 0$, ((2.17)–(2.19)) directly degenerate to those of a two-layer fluid system with a free surface. For this limiting case, these results can also be derived directly from governing equations and boundary conditions and are presented in § A.1.
3. Numerical methods
3.1. A spectral method and DtN operators
The linear model described in ((2.9)–(2.12)) will be numerically solved in spectral space by assuming periodic boundary conditions in horizontal directions on interfacial waves (see figure 1). The spatial Fourier transform in the horizontal plane is defined as
\begin{equation} \mathcal{F}[f](\boldsymbol{k})=\int_{\boldsymbol{x}\in \varOmega}\, {\rm e}^{-\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}f(\boldsymbol{x})\, {\rm d}{\kern0.8pt}\boldsymbol{x}, \end{equation}
where 
$\boldsymbol {x}=(x,y)$ denotes the horizontal variables, 
$\boldsymbol {k}=(k_x,k_y)$ represents the wavevector with 
$k=\sqrt {k_x^2+k_y^2}$ the corresponding wavenumber, 
$\operatorname {\mathrm {i}}=\sqrt {-1}$ is the imaginary unit and 
$\varOmega =[0,L]\times [0,L]$ is the physical domain for integration, namely a region in the horizontal plane with width and length as 
$L$. The set of wavevectors is defined as 
$\boldsymbol {k}\in \varLambda =\mathbb {Z}^2 (2{\rm \pi} /L)$ and 
$\varLambda ^* =\varLambda \setminus {\boldsymbol {0}}$ is defined for the set excluding the zero wavenumber element.
The algorithm consists primarily of three parts: (i) deriving the expression of the DtN operators from the elliptic boundary value problems, (ii) constructing the time-evolution equations, (iii) calculating the topographic coefficients (when the bottom has a spatially varying form). In the subsequent analyses, the particular case with a uniform depth in the lower layer is considered in the first place, where the DtN operators will be defined, and how they work to simulate the evolution of interfacial waves in a three-layer system will be presented. After gaining an overall understanding of the calculation process, we then introduce the bottom topography and elaborate on the more complicated DtN operators.
3.1.1. Elliptic boundary value problems and time-evolution equations
The whole linear system can be divided into three elliptic boundary value problems and one time-evolution problem, and their connection can be achieved by introducing the DtN operators. The three elliptic boundary value problems associated with three-fluid layers are, for the bottom layer,
\begin{equation} \left.\begin{array}{c} \mu^2(\phi_{1xx}+\phi_{1yy})+\phi_{1zz}=0 \quad \text{for } b< z<1,\\ \phi_{1z}-\mu^2(\phi_{1x}b_x+\phi_{1y}b_y)=0 \quad \text{on } z=b,\\ \phi_1=q_1 \quad \text{on } z=1, \end{array} \right\} \end{equation}for the middle layer,
\begin{equation} \left.\begin{array}{c} \mu^2(\phi_{2xx}+\phi_{2yy})+\phi_{2zz}=0 \quad \text{for } 1< z<1+\varTheta_1,\\ \phi_2=q_{22} \quad \text{on } z=1+\varTheta_1,\\ \phi_2=q_{12} \quad \text{on } z=1, \end{array} \right\} \end{equation}and for the upper layer,
\begin{equation} \left.\begin{array}{c} \mu^2(\phi_{3xx}+\phi_{3yy})+\phi_{3zz}=0 \quad \text{for } 1+\varTheta_1< z<1+\varTheta_1+\varTheta_2,\\ \phi_3=q_3 \quad \text{on } z=1+\varTheta_1,\\ \phi_{3z}=0 \quad \text{on } z=1+\varTheta_1+\varTheta_2, \end{array} \right\} \end{equation}
where we introduce the interfacial velocity potentials 
$q_1(\boldsymbol {x},t)=\phi _1(\boldsymbol {x},1+\epsilon \eta _1(\boldsymbol {x},t),t)$, 
$q_{12}(\boldsymbol {x},t)=\phi _2(\boldsymbol {x},1+\epsilon \eta _1(\boldsymbol {x},t),t)$, 
$q_{22}(\boldsymbol {x},t)=\phi _2(\boldsymbol {x},1+\varTheta _1+ \epsilon \eta _2(\boldsymbol {x},t),t)$ and 
$q_3(\boldsymbol {x},t)=\phi _3(\boldsymbol {x},1+\varTheta _1+ \epsilon \eta _2(\boldsymbol {x},t),t)$ as new variables. We show that 
$q_{12}$ and 
$q_{22}$ can be expressed with 
$q_1$ and 
$q_3$ due to the kinematic boundary conditions on two interfaces. Furthermore, we remark that the boundary value problem corresponding to the middle layer has two Dirichlet conditions rather than one Dirichlet and one Neumann in the other two layers. The matchings of solutions on two interfaces can be described as
\begin{equation} \left.\begin{array}{c} \phi_{1z}(\boldsymbol{x},z=1,t)=\phi_{2z}(\boldsymbol{x},z=1,t),\\ \phi_{2z}(\boldsymbol{x},z=1+\varTheta_1,t)=\phi_{3z} (\boldsymbol{x},z=1+\varTheta_1,t), \end{array} \right\} \end{equation}
leading to the elimination of 
$q_{12}$ and 
$q_{22}$ with 
$q_{12}=\mathcal {H}_1[q_1, q_3]$ and 
$q_{22}=\mathcal {H}_2[q_1, q_3]$, where 
$\mathcal {H}_1$ and 
$\mathcal {H}_2$ are pseudo-differential operators defined later on.
The kinematic and dynamic boundary conditions on the interfaces give the time evolutions of interfacial variables, more specifically,
\begin{equation} \left.\begin{array}{c} \eta_{1t}=\dfrac{1}{\mu^2}\mathcal{G}_1[q_1], \quad \eta_{2t}=\dfrac{1}{\mu^2}\mathcal{G}_3[q_3], \\ q_{1t}-R_1q_{12t}=(R_1-1)\eta_1,\quad q_{22t}-R_2q_{3t}=(R_2-1)\eta_2, \end{array} \right\} \end{equation}
where 
$\mathcal {G}_1$ and 
$\mathcal {G}_3$ are the DtN operators defined as 
$\mathcal {G}_1[q_1]=\phi _{1z}(\boldsymbol {x},z=1,t)=\phi _{2z}(\boldsymbol {x},z=1,t)$ and 
$\mathcal {G}_3[q_3]=\phi _{3z}(\boldsymbol {x},z=1+\varTheta _1,t)$, respectively. Once the DtN operators can be appropriately computed, the numerical procedure is summarized as follows: (i) simulation starts with initial conditions for interfacial velocity potentials (
$q_1$ and 
$q_3$) and wave profiles (
$\eta _1$ and 
$\eta _2$); (ii) elliptic boundary problems ((3.2)–(3.4)) are solved, and expressions of the DtN operators are obtained; (iii) time stepping based on (3.6) is performed with the explicit forth-order Runge–Kutta scheme; (iv) simulation repeats the same three steps precedent in the next time step. Apparently, for a three-dimensional problem, the calculation for wave evolution has been reduced to a two-dimensional problem, benefiting from the introduction of the DtN operators. In the subsequent sections we show that this numerical method can be easily implemented and performed in spectral space.
3.1.2. The flat bottom case
In this section the case of trivial bottom topography (i.e. flat bottom) in the lower layer is considered. The formulations of the DtN operators in this situation are given by conducting mathematical derivations based on the spectra of the velocity potentials and wave displacements. Considering 
$b(\boldsymbol {x})=0$, all the elliptic boundary value problems become linear, allowing the application of the Fourier transform. The spectral variables are assigned with a hat symbol. Solving the Laplace equations in the spectral space along with boundary conditions, one obtains, in the lower layer,
\begin{equation} \left.\begin{array}{c} \displaystyle \phi_1(\boldsymbol{x},z)=\sum_{\boldsymbol{k}\in\varLambda}\, {\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\widehat{q_1} (\boldsymbol{k})\dfrac{\cosh(\mu kz)}{\cosh(\mu k)},\\ \displaystyle \phi_{1z}(\boldsymbol{x},z=1)=\sum_{\boldsymbol{k}\in\varLambda}\, {\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\widehat{q_1}(\boldsymbol{k})\mu k\tanh(\mu k), \end{array} \right\} \end{equation}
where 
$\widehat {q_1}(\boldsymbol {k})$ is the spectrum, namely the spatial Fourier transform of 
$q_1(\boldsymbol {x},t)$ at specific time 
$t$. Similarly, in the upper layer the velocity potential and its normal derivative on the upper interface can be expressed as
\begin{equation} \left.\begin{array}{c} \displaystyle \phi_3(\boldsymbol{x},z)=\sum_{\boldsymbol{k}\in\varLambda}\, {\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\widehat{q_3}(\boldsymbol{k})\dfrac{\cosh(\mu k(z-1-\varTheta_1-\varTheta_2))}{\cosh(\mu k\varTheta_2)},\\ \displaystyle \phi_{3z}(\boldsymbol{x},z=1+\varTheta_1)=\sum_{\boldsymbol{k}\in\varLambda}\, {\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\widehat{q_3}(\boldsymbol{k})[-\mu k\tanh(\mu k\varTheta_2)]. \end{array} \right\} \end{equation}As noted in the previous section, the middle-layer fluid is bounded by two interfaces, which leads to two Dirichlet boundary conditions. The solution form of the velocity potential is slightly different, namely
\begin{equation} \left.\begin{array}{c} \displaystyle \phi_2(\boldsymbol{x},z)=\sum_{\boldsymbol{k}\in\varLambda^*}\, {\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\left[\widehat{q_{22}}(\boldsymbol{k})\dfrac{\sinh\left(\mu k(z-1)\right)}{\sinh(\mu k\varTheta_1)}+\widehat{q_{12}}(\boldsymbol{k})\dfrac{\sinh\left(\mu k(1+\varTheta_1-z)\right)}{\sinh(\mu k\varTheta_1)}\right],\\ \displaystyle \phi_{2z}(\boldsymbol{x},z=1)=\sum_{\boldsymbol{k}\in\varLambda^*}\, {\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\left[\widehat{q_{22}}(\boldsymbol{k})\dfrac{\mu k}{\sinh(\mu k\varTheta_1)}-\widehat{q_{12}}(\boldsymbol{k})\dfrac{\mu k}{\tanh(\mu k\varTheta_1)}\right],\\ \displaystyle \phi_{2z}(\boldsymbol{x},z=1+\varTheta)=\sum_{\boldsymbol{k}\in\varLambda^*}\, {\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\left[\widehat{q_{22}}(\boldsymbol{k})\dfrac{\mu k}{\tanh(\mu k\varTheta_1)}-\widehat{q_{12}}(\boldsymbol{k})\dfrac{\mu k}{\sinh(\mu k \varTheta_1)}\right]. \end{array} \right\} \end{equation}
We remark that the spectrum of 
$\phi _2$ at the zero wavevector should vanish. In fact, 
$\phi _2$ can be recast to a new form, 
$\phi _2(x,y,z)=\bar {\phi }_2(z)+\phi _2^*(x,y,z)$, where the perturbation part with zero mean on a horizontal plane at any height, 
$\phi _2^*(x,y,z)$, and the background uniform part, 
$\bar {\phi }_2(z)$, can be separated. With the model being linear, the superposition rule holds, and from the Laplace equation, the horizontally uniform part writes
\begin{equation} \bar{\phi}_2(z)=Az+B, \quad \bar{\phi}_{2z}(z)=A. \end{equation}
The continuity of the vertical velocity component on the lower interface (see (3.5)) requires 
$A=\mathcal {F}[\phi _{1z}(\boldsymbol {x},z=1)](\boldsymbol {0})=0$. We can choose 
$B=0$ without loss of generality and take 
$\widehat {q_{22}}(\boldsymbol {0})=\widehat {q_{12}}(\boldsymbol {0})=\widehat {\phi _2}(\boldsymbol {0})=0$. On the other hand, it is reasonable to apply 
$\widehat {q_1}(\boldsymbol {0})=\widehat {q_3}(\boldsymbol {0})=0=\widehat {\eta _1}(\boldsymbol {0})=\widehat {\eta _2}(\boldsymbol {0})$, and in numerical simulations the time evolution of spectra on wavevector 
$\boldsymbol {k}=\boldsymbol {0}$ can be ignored.
With the aid of the orthogonality of the Fourier primary functions, for wavevector 
$\boldsymbol {k}\in \varLambda ^*$, the matching conditions in (3.5) lead to
\begin{equation} \left.\begin{array}{c}
\widehat{q_{12}}(\boldsymbol{k})={-}\widehat{q_1}(\boldsymbol{k})\tanh(\mu
k)\tanh(\mu
k\varTheta_1)+\widehat{q_{22}}(\boldsymbol{k})\dfrac{1}{\cosh(\mu
k\varTheta_1)},\\
\widehat{q_{12}}(\boldsymbol{k})=\widehat{q_3}(\boldsymbol{k})\sinh(\mu
k\varTheta_1)\tanh(\mu
k\varTheta_2)+\widehat{q_{22}}(\boldsymbol{k})\cosh(\mu k
\varTheta_1). \end{array} \right\}
\end{equation}
It directly follows that 
$\widehat {q_{12}}$ and 
$\widehat {q_{22}}$ can be solved from the above linear algebra system and expressed with 
$\widehat {q_1}$ and 
$\widehat {q_3}$.
The DtN operators and the time-evolution equations in the spectral space can then be expressed as
\begin{equation} \left.\begin{array}{c} \widehat{\mathcal{G}_1[q_1]}(\boldsymbol{k})=\mathcal{F}[\phi_{1z} (\boldsymbol{x},z=1)](\boldsymbol{k})=\widehat{q_1} (\boldsymbol{k})\mu k\tanh(\mu k),\\ \widehat{\mathcal{G}_3[q_3]}(\boldsymbol{k})=\mathcal{F} [\phi_{3z}(\boldsymbol{x},z=1+\varTheta_1)](\boldsymbol{k}) ={-}\widehat{q_3}(\boldsymbol{k})\mu k\tanh(\mu k\varTheta_2),\\ \widehat{\eta_1}_t(\boldsymbol{k})=\dfrac{1}{\mu^2}\widehat{\mathcal{G}_1[q_1]} (\boldsymbol{k}),\quad \widehat{q_1}_t(\boldsymbol{k})-R_1\widehat{q_{12}}_t (\boldsymbol{k})=(R_1-1)\widehat{\eta_1}(\boldsymbol{k}),\\ \widehat{\eta_2}_t(\boldsymbol{k})=\dfrac{1}{\mu^2}\widehat{\mathcal{G}_3[q_3]} (\boldsymbol{k}),\quad\widehat{q_{22}}_t(\boldsymbol{k})-R_2\widehat{q_3}_t (\boldsymbol{k})=(R_2-1)\widehat{\eta_2}(\boldsymbol{k}). \end{array} \right\} \end{equation}Therefore, numerical calculations of time stepping in (3.12) can be effectuated easily without matrix inversion operations.
3.2. Topographic effects
Next, we consider a more general case when the bottom topography is non-trivial. Inspired by Andrade & Nachbin (Reference Andrade and Nachbin2018), the linear three-dimensional formulation is used to capture the effects of large and rapidly varying bottom topography. The variable topography is involved in the elliptic boundary value problem of the bottom layer. Consequently, only in the solution to (3.2) do we need to introduce the topographic coefficient 
$X_{\boldsymbol {k}}$ in the spectral space. It is worth pointing out that Craig, Guyenne & Selum (Reference Craig, Guyenne and Selum2015) have also noted this analogy for the topographic contribution to the DtN operator between the one-layer and two-layer settings. We will briefly describe the non-local formulation and Galerkin approximation used for solving topographic coefficients since they are similar to the derivations of the single-layer case presented by Andrade & Nachbin (Reference Andrade and Nachbin2018).
Following Milewski (Reference Milewski1998) or Andrade & Nachbin (Reference Andrade and Nachbin2018), the velocity potential on the bottom layer can be expressed as
\begin{equation} \phi_1(\boldsymbol{x},z)=\sum_{\boldsymbol{k}\in\varLambda^*}\, {\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}} \left[\widehat{q_1}(\boldsymbol{k})\frac{\cosh(\mu kz)}{\cosh(\mu k)}+X_{\boldsymbol{k}}\frac{\sinh(\mu k(z-1))}{k\cosh^2(\mu k)}\right]. \end{equation}
Substituting (3.13) into the impermeability boundary condition at the bottom, the relation between topographic coefficients and the Dirichlet data 
$q_1$ is established, namely,
\begin{equation} \sum_{\boldsymbol{k}\in\varLambda^*}\widehat{q_1}(\boldsymbol{k}) \boldsymbol{\nabla} \boldsymbol{\cdot} \left[{\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\frac{-\sinh(\mu kb)}{\cosh(\mu k)}\frac{\boldsymbol{k}}{k}\right]= \sum_{\boldsymbol{k}\in\varLambda^*}X_{\boldsymbol{k}}\boldsymbol{\nabla} \boldsymbol{\cdot} \left[{\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\frac{\cosh(\mu k(b-1))}{\cosh^2(\mu k)}\frac{\boldsymbol{k}}{k^2}\right]. \end{equation}
The wavenumber truncation is achieved by the Galerkin approximation to reduce the computational costs and increase the stability of the algorithm. For a numerical simulation with 
$N^2$ grid points, 
$N^2$ modes exist in the spectral space. For an arbitrary real function 
$f$,
\begin{equation} \hat{f}(\boldsymbol{k})=\overline{\hat{f}(-\boldsymbol{k})} \end{equation}
holds, where the overbar symbol stands for complex conjugation; thus, only half of the modes need to be calculated. The maximum wavenumber is 
$N\sqrt {2}{\rm \pi} /L$. Under the Galerkin approximation, we only calculate topographic coefficients on wavevectors whose modulus is less than a critical value 
$k_m$, and the number of these wavevectors is 
$M \ll N^2/2$. Coefficients on wavevectors with larger modulus are set to 0, assuming that the interaction between these components and interfacial waves remains negligible. Then, the corresponding set of wavevectors is defined as 
$\varLambda _M^*=\{\boldsymbol {k} \in \varLambda ^* | k< k_m\}$.
We define two linear operators that map a sequence of complex coefficients to a function of 
$\boldsymbol {x}$ as
\begin{equation} \left.\begin{array}{c} \displaystyle \mathcal{A}[f_{\boldsymbol{k}}]= \sum_{\boldsymbol{k}\in\varLambda^*}-f_{\boldsymbol{k}} \boldsymbol{\nabla} \boldsymbol{\cdot} \left[{\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\dfrac{\sinh(\mu kb)}{\cosh(\mu k)}\dfrac{\boldsymbol{k}}{k}\right], \\ \displaystyle \mathcal{B}_M[f_{\boldsymbol{k}}]= \sum_{\boldsymbol{k}\in\varLambda_M^*}f_{\boldsymbol{k}} \boldsymbol{\nabla} \boldsymbol{\cdot} \left[{\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\dfrac{\cosh(\mu k(b-1))}{\cosh^2(\mu k)} \dfrac{\boldsymbol{k}}{k^2}\right]. \end{array} \right\} \end{equation}
Here 
$\mathcal {A}[\widehat {q_1}(\boldsymbol {k})]$ is a smooth function, denoted as 
$Y(\boldsymbol {x})$, provided 
$\widehat {q_1}(\boldsymbol {k})$ is known. The Galerkin approximation requires that the residual 
$R_M=\mathcal {B}_M[X_{\boldsymbol {k}}]-Y$ satisfies
\begin{equation} (R_M,{\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{l}\boldsymbol{\cdot} \boldsymbol{x}})=0 \quad\text{for } \forall\boldsymbol{l}\in\varLambda_M^*, \end{equation}which is equivalent to
\begin{align} &\sum_{\boldsymbol{k}\in\varLambda_M^*}X_{\boldsymbol{k}}\frac{\operatorname{\mathrm{i}}(\boldsymbol{l}\boldsymbol{\cdot} \boldsymbol{k})}{k^2}\mathcal{F}\left[{\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\frac{\cosh(\mu k(b-1))}{\cosh^2(\mu k)}\right](\boldsymbol{l})\nonumber\\ &\quad =\sum_{\boldsymbol{k}\in\varLambda^*}\widehat{q_1}(\boldsymbol{k})\frac{-\operatorname{\mathrm{i}}(\boldsymbol{l}\boldsymbol{\cdot} \boldsymbol{k})}{k}\mathcal{F}\left[{\rm e}^{\operatorname{\mathrm{i}}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\frac{\sinh(\mu kb)}{\cosh(\mu k)}\right](\boldsymbol{l}). \end{align}
Equation (3.18) can be transformed into a linear system. The matrix system is constructed by extracting the real and imaginary parts of the complex system. Then, the real matrix will be inverted. Due to the large dimension of the matrix, the computational cost can easily become extremely high and force restrictions on the accuracy of the simulation. The topographic coefficients and their dependence on 
$\widehat {q_1}(\boldsymbol {k})$, i.e. 
$\textrm {d}X_{\boldsymbol {k}_i}/\textrm {d}\widehat {q_1}(\boldsymbol {k}_j)$, can be solved by inverting the system. We remark that the matrix inversion is effectuated with the LU decomposition method, and the dimension of the real matrix constructed and inverted is 
$2(M-1)\times 2(M-1)$. Each spectral component is decomposed into real and imaginary parts, and the zero wavevector is excluded. The dimension of the Jacobian matrix 
$J_{ij}=\textrm {d}X_{\boldsymbol {k}_i}/\textrm {d}\widehat {q_1}(\boldsymbol {k}_j)$ is 
$2(M-1)\times 2(N^2/2-1)$. The validation and convergence study of the Galerkin approximation can be found in Andrade & Nachbin (Reference Andrade and Nachbin2018), as well as the method of choosing the parameter 
$M$. Typically, 
$M$ is selected by considering the maximal amplitude of bottom topography and particle trajectories beneath the free surface in the single-layer fluid case. However, the dimension of the matrix can still easily become very large, so the matrices cannot be stored for the simulation of a three-dimensional system. Therefore, a smaller 
$M$ is used in three-dimensional cases. We also remark that all modes in the spectral space are solved without the Galerkin approximation for two-dimensional cases.
It is worth noting that solutions of the velocity potentials in the boundary value problems for the middle- and upper-layer fluids remain unchanged. Using the expression of 
$\phi _1$ under the non-local formulation in (3.13) and matching solutions on two interfaces (see (3.5)), one obtains
\begin{equation} \left.\begin{array}{c}
\displaystyle
\widehat{q_{12}}(\boldsymbol{k})={-}\dfrac{\tanh(\mu
k)}{\tanh(\mu
k\varTheta_1)}\widehat{q_1}(\boldsymbol{k})-\dfrac{1}{k\cosh^2(\mu
k)\tanh(\mu
k\varTheta_1)}X_{\boldsymbol{k}}-\dfrac{\tanh(\mu
k\varTheta_2)}{\sinh(\mu
k\varTheta_1)}\widehat{q_3}(\boldsymbol{k}),\\
\displaystyle
\widehat{q_{22}}(\boldsymbol{k})={-}\dfrac{\tanh(\mu
k)}{\sinh(\mu
k\varTheta_1)}\widehat{q_1}(\boldsymbol{k})-\dfrac{1}{k\cosh^2(\mu
k)\sinh(\mu
k\varTheta_1)}X_{\boldsymbol{k}}-\dfrac{\tanh(\mu
k\varTheta_2)}{\tanh(\mu
k\varTheta_1)}\widehat{q_3}(\boldsymbol{k}). \end{array}
\right\}
\end{equation}
Substituting the above expressions into (3.6), the time-evolution equations for 
$\boldsymbol {k}\in \varLambda ^*$ read
\begin{equation} \left.\begin{array}{c} \displaystyle\widehat{\eta_1}_{t}(\boldsymbol{k})= \dfrac{1}{\mu^2}\widehat{\mathcal{G}_1[q_1]}(\boldsymbol{k})= \dfrac{1}{\mu}\left[\widehat{q_1}(\boldsymbol{k})k\tanh(\mu k)+X_{\boldsymbol{k}}\dfrac{1}{\cosh^2(\mu k)}\right],\\ \displaystyle\widehat{\eta_2}_{t}(\boldsymbol{k})= \dfrac{1}{\mu^2}\widehat{\mathcal{G}_3[q_3]}(\boldsymbol{k})={-} \dfrac{1}{\mu}\widehat{q_3}(\boldsymbol{k})k\tanh(\mu k\varTheta_2), \end{array} \right\} \end{equation}and
\begin{equation} \left.\begin{array}{c@{}}
\begin{aligned} \widehat{q_3}_{t}(\boldsymbol{k}) & =
\dfrac{\cosh(\mu k\varTheta_1)}{R_1R_2}\left[1+R_1\tanh(\mu
k)\tanh(\mu
k\varTheta_1)\right]\widehat{q_1}_{t}(\boldsymbol{k})\\ &
\quad +\dfrac{\sinh(\mu k\varTheta_1)}{R_2k\cosh^2(\mu
k)}X_{\boldsymbol{k}t}+\dfrac{1-R_1}{R_1R_2}\cosh(\mu
k\varTheta_1)\widehat{\eta_1}(\boldsymbol{k})+
\dfrac{1-R_2}{R_2}\widehat{\eta_2}(\boldsymbol{k}),
\end{aligned}\\[-5pt]\\ \begin{aligned}
\widehat{q_1}_{t}(\boldsymbol{k}) & =
\dfrac{\left(-\dfrac{R_1\tanh(\mu k\varTheta_2)\tanh(\mu
k\varTheta_1)}{R_2k\cosh^2(\mu k)}-\dfrac{R_1}{k\cosh^2(\mu
k)}\right)X_{\boldsymbol{k}t}+\dfrac{R_1(R_2-1)\tanh(\mu
k\varTheta_2)}{R_2\cosh(\mu
k\varTheta_1)}\widehat{\eta_2}(\boldsymbol{k})}{\tanh(\mu
k\varTheta_1)+R_1\tanh(\mu k)+\dfrac{\tanh(\mu
k\varTheta_2)}{R_2}\left[1+R_1\tanh(\mu k)\tanh(\mu
k\varTheta_1)\right]}\\ & \quad
+\dfrac{\left[(R_1-1)\tanh(\mu
k\varTheta_1)+\dfrac{R_1-1}{R_2}\tanh(\mu
k\varTheta_2)\right]\widehat{\eta_1}(\boldsymbol{k})}{\tanh(\mu
k\varTheta_1)+R_1\tanh(\mu k)+\dfrac{\tanh(\mu
k\varTheta_2)}{R_2}\left[1+R_1\tanh(\mu k)\tanh(\mu
k\varTheta_1)\right]}. \end{aligned} \end{array} \right\}
\end{equation}Finally, the time derivative of the topographic coefficient writes
\begin{equation} X_{\boldsymbol{k}t}=\sum_{\boldsymbol{k}_j\in\varLambda_M^*}\frac{\partial X_{\boldsymbol{k}}}{\partial\widehat{q_1}(\boldsymbol{k}_j)} \widehat{q_1}_t(\boldsymbol{k}_j). \end{equation}3.3. Lateral boundary conditions and initialization
In addition to the periodic conditions given implicitly by the spectral method, the relaxation zone method is applied for wave making and absorption in specific areas. This method was successfully used by Guyenne & Nicholls (Reference Guyenne and Nicholls2008) for free-surface water-wave problems and has proven efficient and easy to implement. At each time step, the updated spectra will be transformed backward to the physical space, and the physical variables will then be relaxed towards a specified analytical function over a limited region.
The relaxation is achieved by defining a relaxation coefficient function and a desired solution. In the three-layer system the desired solution can be designated as 
$(\widetilde {q_1}$, 
$\widetilde {q_3}$, 
$\widetilde {\eta _1}$, 
$\widetilde {\eta _2})$. An example of a relaxation coefficient as a function of 
$x$ is
\begin{equation} c_r(x)=\frac{1}{2}+\frac{1}{2}\tanh\left[2{\rm \pi} \left(\frac{x}{l}-\frac{p}{l}\right)\right], \end{equation}
where 
$l$ is the characteristic length and 
$p$ is the central position of the relaxation zone. In more general cases, the physical variables are controlled by
\begin{equation} \begin{pmatrix} q_1(\boldsymbol{x},t)\\ q_3(\boldsymbol{x},t) \\ \eta_1(\boldsymbol{x},t) \\ \eta_2(\boldsymbol{x},t) \end{pmatrix}=c_r(\boldsymbol{x})\begin{pmatrix} q_1(\boldsymbol{x},t)\\ q_3(\boldsymbol{x},t) \\ \eta_1(\boldsymbol{x},t) \\ \eta_2(\boldsymbol{x},t) \end{pmatrix}+[1-c_r(\boldsymbol{x})]\begin{pmatrix} \widetilde{q_1}(\boldsymbol{x},t)\\ \widetilde{q_3}(\boldsymbol{x},t) \\ \widetilde{\eta_1}(\boldsymbol{x},t) \\ \widetilde{\eta_2}(\boldsymbol{x},t) \end{pmatrix}. \end{equation}Figure 4 illustrates the generation of incident monochromatic waves in two-dimensional fluid systems. Mode 1 and mode 2 waves are respectively shown in figures 4(a) and 4(b). This method has proved reliable: wave profiles outside the relaxation zone can maintain the desired amplitude with negligible error, and in the absorbing zone, waves are efficiently dissipated without obvious reflection after some time of evolution. In the next section, incident waves will be generated with wave-making zones in both two- and three-dimensional cases.
Figure 4. Two-dimensional monochromatic interfacial waves generated with relaxation zone: waves on the lower interface (solid curves), waves on the upper interface (dotted curves) and theoretical wave amplitudes (horizontal lines). The system parameters are 
$R_1=R_2=0.6$, 
$\varTheta _1=\varTheta _2=1$ and 
$\mu =0.1$. The wave-absorbing zones are 
$x<5$ and 
$x>34$, and the wave-making zone is 
$5< x<15$.
4. Results
4.1. Two-dimensional internal wave problems
In this part, numerical results for two-dimensional cases are presented. The grid parameters are chosen as 
$\varDelta _x=0.02$ and 
$\varDelta _t=0.01$, satisfying the Courant–Friedrichs–Lewy condition, i.e. 
$\varDelta _x/\varDelta _t>c_\phi$. The linear system is fully solved with 
$M=N$ as the computational cost is affordable for the two-dimensional problem.
4.1.1. Linear wave shoaling
First, we calculate the two-dimensional wave shoaling problem and compare the results with the linear theoretical predictions to validate our numerical procedure. The topography is chosen to vary gradually from a lower level to a higher level, which leads to the gradual reduction of the depth of the bottom layer. The variable-depth profile is similar to the example shown in figure 1, defined by
\begin{equation} b(x)=1-b_0\left[1+\tanh{\left(\frac{x-x_0}{2}\right)}\right]. \end{equation}
We remark that 
$2b_0$ is the amplitude of the bottom undulation and 
$x_0$ is the topographic centre, chosen as 
$L/2$ with 
$L$ the length of the computational domain. For a monochromatic incident wave past the variable bathymetry, the frequency remains invariant and the total energy of the system is conserved, namely,
\begin{equation} \mathcal{E}=[a_2^2(\rho_2-\rho_3)+a_1^2(\rho_1-\rho_2)]c_g= \text{const.}, \end{equation}
where 
$c_g$ stands for the group velocity. On the contrary, the change in depth of the fluid layer can lead to a change in the dispersion relation; equivalently, the wavenumber and 
$a_2/a_1$ will evolve. The group velocity can be calculated from the dispersion relation and then the wave elevation can be theoretically evaluated using (4.2).
The amplitude of the bottom topography is chosen to be substantial and, for cases presented in this section, variation in amplitude is fixed as 
$2b_0=0.75$. Wave profiles that reach the steady state are shown in figures 5–7 for three-layer systems with 
$\varTheta _1=\varTheta _2=1$ and different density ratios. The theoretical predictions of wave elevation on the lower and upper interfaces are drawn with a dashed line and a dash-dotted line, respectively, indicating that the numerical results are in good coherence with the linear theory. In figure 5(b,c) the grid convergence and accuracy of the mode-1 case in figure 5(a) are presented with a quantitative comparison between theoretical prediction and statistical numerical results. The time evolution of the Hamiltonian for this three-layer system is investigated, which increases initially because of wave making and finally arrives at a constant level as it should be an intrinsic conserved quantity (see figure 5b). The relative difference in downstream wave amplitude between the theoretical prediction and averaged numerical results is less than 2 %, and it remains almost the same level with growing grid numbers. As a result, our numerical method is validated. For both modes, the monochromatic waves propagate from left to right with a subtle increase in wavenumber. For a mode 1 wave (see figure 5a), the wave amplitude on the lower interface decreases and the one on the upper interface increases; on the contrary, for a mode 2 wave, the wave amplitude on the lower interface rises and the one on the upper interface falls. This phenomenon differs from the single-layer case calculated by Guyenne & Nicholls (Reference Guyenne and Nicholls2008). In a single-layer system the total energy depends on the wave amplitude and the group velocity with a constant fluid density. The reduction of fluid depth leads to the increase of wavenumber and decrease of group velocity, giving rise to the increase of wave amplitude.
Figure 5. Mode 1 and mode 2 linear waves shoaling over slowly varying bottom topography. The wave-making zone is located in 
$x<12$, and the system parameters are 
$\varTheta _1=\varTheta _2=1$, 
$R_1=0.6$, 
$R_2=0.4$ and 
$\mu =0.1$. (a) Wave profiles: waves on the lower interface (solid lines), waves on the upper interface (dotted lines) and theoretical predictions of wave envelope (dashed and dash-dotted lines). (b) Evolution of Hamiltonian computed with various numbers of collocation points: 
$N=2000$ (black solid line), 
$N=1600$ (blue circles), 
$N=1200$ (red squares) and 
$N=800$ (pink triangles). (c) Grid convergence test: a comparison between the theoretical prediction of downstream wave amplitude and averaged results in numerical simulation.
Figure 6. Mode 1 and mode 2 linear waves shoaling over slowly varying bottom topography. The wave-making zone is located in 
$x<12$. The system parameters are 
$\varTheta _1=\varTheta _2=1$, 
$\mu =0.1$ and various density ratios: (a) 
$R_1=R_2=0.6$ and (b) 
$R_1=0.6$, 
$R_2=0.85$. Solid lines, waves on the lower interface; dotted lines, waves on the upper interface; dashed and dash-dotted lines, theoretical predictions of wave envelope.
Figure 7. A mode 1 linear wave shoaling over slowly varying bottom topography. The wave-making zone is located in 
$x<12$, and the system parameters are chosen as 
$R_1=0.6$, 
$R_2=0.8$, 
$\varTheta _1=\varTheta _2=1$ and 
$\mu =0.1$. Solid line, wave on the lower interface; dotted line, wave on the upper interface; dashed and dash-dotted lines, theoretical predictions of wave envelope.
It can be observed from figures 5 and 6 that the evolution of wave amplitude is more complicated when it comes to a three-layer system. The density ratios influence the change in wave amplitude. Compared with the case with 
$R_2=0.6$ in figure 6(a), figure 5(a) shows that for mode 1 waves when the upper-layer fluid becomes lighter (
$R_2=0.4$), the variation of wave amplitude on the upper interface becomes less evident, while the variation on the lower interface shows no apparent change. The amplitude evolution of a mode 2 wave in figure 5(a) is akin to the case shown in figure 6(a) but with a smaller relative change. As the density of the upper fluid increases, the tendency of wave amplitude evolution on the lower interface changes after 
$R_2$ reaches a critical value of around 
$0.8$. In figure 6(b), with 
$R_2=0.85$, the amplitude of the mode 1 internal wave on the lower interface rises to another level while climbing the bottom topography instead of decreasing in the other two cases with 
$R_2=0.4$ and 
$R_2=0.6$. Conversely, no evident change in the behaviour of the mode 2 wave is observed. In figure 7, near the critical value of 
$R_2$ for mode 1 waves, the wave amplitude on the lower interface rises first and then falls to almost the same level as the incident waves.
The influence of depth ratios is investigated by the theoretical prediction of amplitude evolution curves demonstrated in figure 8. The variation tendency does not seem to change in a few limited benchmark tests. It can be observed that with a thick middle layer, the difference between wave amplitudes on the two interfaces increases for mode 1 and decreases for mode 2. This phenomenon also appears when the upper layer becomes thicker. We noticed that in the thick middle-layer configuration, the wave amplitude on the lower interface reaches a level close to zero for mode 1. For mode 2, the evanescent waves are present on the upper interface. Besides, the change in thickness of the upper layer exerts a smaller influence on waves on the lower interface.
Figure 8. Theoretical predictions of wave amplitude for 
$R_1=R_2=0.6$, 
$\mu =0.1$ and different values of depth ratio: 
$\varTheta _1=\varTheta _2=1$ (solid line); 
$\varTheta _1=4$, 
$\varTheta _2=1$ (dashed line); 
$\varTheta _1=1$, 
$\varTheta _2=4$ (dotted line). The wave amplitude on the lower interface is presented with a single line, and on the upper interface it is presented with a marked line.
4.1.2. Gaussian wavepacket past a significant obstacle
We replace the monochromatic incident wave with a wavepacket of Gaussian envelope. The objective is to observe the excitation of internal waves of mode 2 when a mode 1 wavepacket passes over a local bottom topography with a significant height (or in a reverse way). In the case of a flat bottom, the main part of the wavepacket should maintain its original mode. Nevertheless, we show that narrowing the obstacle while keeping its height makes the mode excitation possible in a system with suitable parameters.
The initialization of the incident wave is achieved by filtering a monochromatic wave with a Gaussian window. We take the filter function as
\begin{equation} f_g(x)={\rm e}^{-{(x-x_c)^2}/{2\sigma^2}}, \end{equation}
where 
$x_c$ is the centre of the filtered wavepacket, and the maximal amplitude of the filter is 1 such that the original amplitude of the monochromatic wave is preserved as the maximum of the wavepacket. The bottom obstacle features a half-period cosine-type profile. Mathematically, 
$b(x)$ can be defined as
\begin{equation} b(x)=\left\{ \begin{array}{@{}ll} b_0\cos^2\left(\dfrac{{\rm \pi}(x-x_0)}{L_0}\right), & |x-x_0|<\dfrac{L_0}{2},\\ 0, & |x-x_0|>\dfrac{L_0}{2}, \end{array}\right. \end{equation}
where 
$b_0$ is the height of the bottom obstacle and 
$L_0$ is the width. The centre position of the bottom obstacle is 
$x_0=18$, while the length of the computational domain is set as 
$L=40$. The height is chosen to be significant with 
$b_0=0.8$ for all the cases presented in the following. We vary the width so that the front of the bottom topography has different slopes.
Figure 9 presents evolutions of mode 1 interfacial waves passing over bottom topography with different widths. Each layer has the same depth in this system and the density ratios are 
$R_1=R_2=0.6$, which means the density difference between adjacent layers is relatively evident. The wavepacket extends to a broader range during propagation due to the dispersive effect. When the wavepacket passes a wide obstacle with 
$L_0=8$, no mode change occurs from beginning to end, as shown in figure 9(a), though the wave amplitude varies during the stage when the wave climbs up the bottom topography. With a narrower obstacle, 
$L_0=4$ in figure 9(b), a mode 2 wave gradually appears while passing over the obstacle. The mode 2 wave is ultimately excited and propagates with a slightly higher wavenumber and a slower phase velocity. Consequently, it is mainly observed behind the mode 1 wave in the trailing edge of the packet. With 
$L_0=2$, the excited mode 2 wave has a larger amplitude in figure 10(a) compared with the case in figure 9(b). Moreover, figure 10(b) shows that a narrow bottom obstacle can also excite a mode 1 wave when a mode 2 wavepacket passes over. In this case, the generated mode-1 component surpasses the mode-2 part due to its lower wavenumber and higher phase velocity. Finally, to show that the mode exchanging caused by a locally confined bottom topography is also possible in more realistic situations, we consider a case with both density ratios close to 1, namely 
$R_1=R_2=0.9$, and the depth of each layer remains the same. In figure 11, we observe an evident excitation and propagation of a mode 2 wave, and the whole wavepacket is dispersed into a broader range.
Figure 9. Evolution of an internal wavepacket with a Gaussian envelope past a locally confined, significant bottom topography. The width of the obstacle is 
$L_0=8$ for (a) and 
$L_0=4$ for (b), and the other parameters are 
$\varTheta _1=\varTheta _2=1$ and 
$R_1=R_2=0.6$. Wave profiles at 
$t=0$, 
$t=25$, 
$t=50$ and bottom topography are presented from top to bottom. The height of the obstacle is 
$b_0=0.8$. Solid line, wave on the lower interface; dotted line, wave on the upper interface.
Figure 10. Evolution of internal wavepackets with a Gaussian envelope past a locally confined, significant bottom topography: (a) mode-1 initial profile and (b) mode-2 initial profile. The height of the obstacle is 
$b_0=0.8$, the width is 
$L_0=1$ and the other system parameters are 
$\varTheta _1=\varTheta _2=1$ and 
$R_1=R_2=0.6$. Wave profiles at 
$t=0$, 
$t=25$, 
$t=50$ and bottom topography are presented from top to bottom. Solid line, wave on the lower interface; dotted line, wave on the upper interface.
Figure 11. Evolution of an internal wavepacket with a Gaussian envelope past a locally confined, significant topography. The height of the obstacle is 
$b_0=0.8$, the width is 
$L_0=1$ and the other system parameters are 
$\varTheta _1=\varTheta _2=1$ and 
$R_1=R_2=0.9$. Wave profiles at 
$t=0$, 
$t=50$, 
$t=100$ and bottom topography are presented from top to bottom. The solid and dotted lines represent, respectively, interfacial waves on the lower and upper interfaces.
4.1.3. Quickly varying topography
We next consider another interesting situation where an internal wavepacket interacts with a quickly varying bottom topography. The internal wavepacket is initialized with a carrier wave, a monochromatic incident wave of wavelength 
$\lambda =1$, filtered by a hyperbolic tangent function. Similarly, the bottom obstacle is described by a sinusoidal wave of wavelength 
$l$ multiplied by a hyperbolic tangent window function. The system parameters are fixed to 
$\varTheta _1=\varTheta _2=1$ and 
$R_1=R_2=0.6$. The fundamental wavelength of interfacial waves is also invariant. We investigate how the wavepacket evolves when passing over obstacles with different basic wavelengths.
Numerical experiments have been carried out for two cases: 
$\lambda \approx l$ and 
$\lambda \gg l$. Figures 12 and 13 present wave profiles at different instants as the wavepacket passes over two kinds of bottom topography: quickly varying topography and plateau-shaped topography. The fundamental wavelength of wavepacket bottom topography varies from 
$l=0.8$ to 
$l=0.25$.
Figure 12. Evolution of mode 1 internal wavepackets with a hyperbolic tangent envelope passing over a quickly varying bottom topography. The fundamental wavelength for the wavepackets is 
$\lambda =1$ and, for the bottom topography, it is (a) 
$l=0.8$ and (b) 
$l=0.4$. Wave profiles at 
$t=1$, 
$t=10$, 
$t=20$ and bottom topography are presented from top to bottom. Solid lines, wave on the lower interface; dotted lines, wave on the upper interface.
Figure 13. Evolution of mode 1 internal wavepackets with a hyperbolic tangent envelope passing over a quickly varying bottom topography. The basic wavelength for the wavepackets is 
$\lambda =1$ and, for the bottom topography, it is (a) 
$l=0.267$ and (b) 
$l=0.25$. Wave profiles at 
$t=1$, 
$t=10$, 
$t=20$ and bottom topography are presented from top to bottom. Solid lines, wave on the lower interface; dotted lines, wave on the upper interface.
Before reaching the undulation part of the bottom, the wave envelope deforms slightly due to the dispersive effect. Over the bottom topography, small-amplitude ripples are excited, ascribed to the perturbation applied to the incident wavepacket. These oscillations are close to standing waves and become more evident when the bottom topography varies more quickly. The typical wavelength of these ripples is approximately 
$1.5l$. For comparison, in the case where the bottom topography does not present oscillations, these small standing waves do not appear. We also remark here that when the wavelength of bottom topography is small, the numerical stability is challenged. The grid number and computational domain cannot increase infinitely, which will lead to the divergence of the standing waves over the bottom topography. No dissipation in the numerical algorithm with the spectral method and the difficulty of solving complex topographic effects by matrix inversion may be the reason for the divergence. When the topography varies quickly, the increase in grid point or change in the computational domain will easily lead to different topographic coefficients. The increase in complexity will evidently influence the accuracy of matrix inversion, which is worth further improvement on the treatments.
For internal wavepackets equipped with mode-1 carrier waves interacting with bottom topography, we focus on reflected waves and the potential excitation of mode 2 waves. During the wavepacket evolution, small ripples propagating to the left appear on the trailing edge in the early stage due to the filtering on the incident wave envelope. Then, a reflected wave is generated when the wavepacket meets the bottom obstacle. It is worth noting that the reflected wave may interact with the excited mode 2 wave as the latter has a slower phase speed. In figure 12(a) the reflected wave has negligible amplitude and no mode 2 wave is excited. With the decrease of the fundamental wavelength of bottom topography, in figure 12(b), both a reflected wave and a mode 2 internal wave can be observed in the left to the bottom topography (
$6< x<10$). The wavelength in the domain 
$6< x<8$, where the mode 1 wave dominates, is 
$\lambda _r\approx 0.8=2l$, and the one in the domain 
$8< x<10$, where the mode 2 wave presents, is 
$\lambda _{m2}\approx 0.68<2l$. The first wavelength corresponds to the Bragg resonance, and the second one is smaller than the resonant reflected wavelength 
$2l$ but greater than the theoretical mode-2 wavelength 
$\lambda _{t,m2}\approx 0.625$ excited with incident mode-1 wavelength 
$\lambda =1$. As a result, in the interval 
$[8,10]$, the wave presents a superposition of two components. The characteristics of wave profiles in figure 13(a) in the same area are akin to the precedent analysis with a smaller 
$l$.
From figures 12(a) to 13(b), the wavelength of the bottom topography decreases from 
$l=0.8$ to 
$0.25$. Mode 2 waves are only excited with a limited range of 
$l$. It is shown that in figures 12(b)–13(a) the amplitude of mode 2 waves decreases, and then with an even smaller 
$l$ in figure 13(b), the mode 2 wave disappears. Furthermore, figure 14 illustrates the homogenization effect of long internal waves passing over a quickly varying bottom obstacle. No mode change exists during the propagation and the wave profiles in figure 14 are similar to the case with smooth bottom topography in the envelope form and to the flat bottom case. It can be observed that the wave passing over quickly varying bottom topography has a phase velocity smaller than that passing over the envelope but faster than the flat bottom case. Therefore, the results imply that the internal waves feel the change in depth in a homogenized way, i.e. the internal waves on both interfaces feel the delicate structures of bottom topography in a mean-field approximation. These results demonstrate that the homogenization effect also takes action for internal waves.
Figure 14. A comparison among wave profiles at 
$t=20$ for wavepackets past over different topographies: a quickly varying bottom topography (solid lines), bottom topography with hyperbolic tangent envelope (dashed line with blue circles) and the flat bottom case (dotted lines with red squares). The quickly varying bottom topography is in the shape of a wavepacket with the basic wavelength 
$l=0.222$.
Another particular phenomenon should also be highlighted in our cases. On the left side, far from the bottom topography (
$2< x<6$), resonance exists on the wavelength 
$\lambda _{r2}\approx 4l$. This is clearly observed in figure 13(a), as well as in figures 12(b) and 14 (solid lines). This resonance is emphasized with a decrease of 
$l$ from 
$0.4$ to 
$0.267$ and then weakened with a further reduction from 
$0.267$ to 
$0.222$. The resonance is excited around 
$t\approx 13$ when the incident wavepacket is entirely over the bottom topography and interacts directly with it. The resonant harmonic then travels to the left, similar to reflected waves. The numerical model assumes that the height of the topography is in the same order as the fluid depth of the bottom layer, which leads to a strong coupling between the wave on the lower interface and the topography. Consequently, resonance among harmonics becomes possible.
4.2. Three-dimensional problems
In this section we delve into the refraction of horizontally two-dimensional interfacial waves in response to variations in bathymetry. Our investigation begins with a detailed examination of interfacial waves passing over a specific submerged circular mound, often termed a Luneberg lens. This unique bathymetric feature offers valuable insights into wave refraction phenomena. Subsequently, we shift our focus to the evolution of interfacial waves in conjunction with free-surface waves within a specific scenario, namely a two-and-a-half-layer system. This investigation takes into account the influence of realistic bathymetry, such as the complicated bathymetric patterns found near the Strait of Gibraltar. We aim to comprehensively understand how bathymetric variations impact wave behaviour in these settings.
4.2.1. Refraction by a submerged mound
A circular submerged mound is placed on the bottom of a three-layer system with finite depth. Each layer has equal thickness with 
$\varTheta _1=\varTheta _2=1$. This particular mound plays a role in bending the wavefront in a specific way, and this effect is analogous to that of the Luneburg lens in optics. In addition to the numerical simulation, the ray theory is applied to obtain the theoretical prediction for propagation directions of the wavefront, namely rays, to show that our numerical method can accurately capture the changes in propagation direction in a three-dimensional system. The bottom topography is given by
\begin{equation} b_m(r)=\left\{ \begin{array}{@{}ll} 1-\dfrac{\alpha^2}{\alpha^2+1-(r/r_0)^2} & \text{if } r< r_0,\\ 0 & \text{if } r\ge r_0, \end{array}\right. \end{equation}
where we consider 
$\alpha =0.85$ and 
$r=r(\boldsymbol {x})=\sqrt {(x-x_c)^2+(y-y_c)^2}$; 
$\boldsymbol {x}_c=(x_c,y_c)$ is the centre of the mound and 
$r_0$ stands for the radius. Here we fix 
$x_c=y_c=5.5$ and 
$r_0=1.4$. The computational domain is a square with dimensions 
$L\times L$ (we choose 
$L=12$ in subsequent simulations), and 240 grid points are used in each horizontal direction. Due to the memory restrictions while solving the inversion of the matrix, we applied the Galerkin approximation with the parameter 
$M=14 400$, indicating that we only calculate the topographic coefficients on the first 
$M$ wavevectors with relatively small modulus.
A mode-1 monochromatic incident wave propagates from left to right with the wavelength 
$\lambda =1$. Snapshots for interfacial wave profiles at 
$t=5$ and 
$t=10$ are presented in figure 15, respectively. Before arriving at the region of terrain undulation, the incident wave keeps its original monochromatic form well. When the wavetrain enters the terrain-affected region, the topographic effect becomes active and the wavefronts start to deform. The fluid depth of the bottom layer decreases due to the convex mound, reducing the group velocity. At 
$t=5$, the wavefronts on both interfaces bend backward to the left, forming a semicircular-like structure.
Figure 15. A mode 1 internal wave passing over a convex mound with the Luneburg lens profile: (a) 
$t=5$ and (b) 
$t=10$.
It is worth noting that for mode 1 internal waves, the wave amplitude on the lower interface in the central area is reduced to a lower level than the original amplitude (see figure 16). On the contrary, the wave amplitude presents a higher value on the upper interface. Wave energy on two interfaces may exchange through the middle layer via the kinematic and dynamic boundary conditions. This means interfacial waves on the lower interface can influence those on the upper interface as they are coupled in time-evolution equations.
Figure 16. A comparison between numerical solutions and theoretically predicted rays for a mode-1 incident wave. In the upper figures the dashed lines stand for rays, and the contour plots are the isolines of the internal wave numerically calculated. Wave profiles along the streamwise direction at 
$y=y_c=5.5$ are shown in the lower figures for 
$t=9.5$ (solid lines), 
$t=10$ (dashed lines), 
$t=10.5$ (dotted lines) and 
$t=11$ (dash-dotted lines). (a) Waves on the lower interface (the amplitude reaches its maximum at 
$x\approx 7.45$); (b) waves on the upper interface (the amplitude reaches its maximum at 
$x\approx 8$).
It can be observed in figure 16 that the potential energy on both interfaces is focused on a point downstream, behind the central position of the mound. At 
$t\approx 10$, the amplitude of the mode 1 wave in the vicinity of the focusing point reaches approximately twice the original amplitude. More downstream behind this point, the wave amplitudes decay.
The ray theory predicts the refraction for monochromatic waves and can be used to verify the refraction pattern of the numerical wave fields. Considering weak dispersion with a small parameter 
$\mu$, the ansatz for the velocity potential writes
\begin{equation} \phi_j(\boldsymbol{x},t)=Am_j(\boldsymbol{x})\,{\rm e}^{{\rm i}(\theta(\boldsymbol{x})-\omega t)}, \quad j=1,2, \end{equation}
where 
$Am_j(\boldsymbol {x})$ is an amplitude function and 
$\theta (\boldsymbol {x})$ is a phase function. Both 
$Am_j(\boldsymbol {x})$ and 
$\theta (\boldsymbol {x})$ vary faster than the depth variation of the lower-layer fluid. The eikonal equation governs the spatial distribution of the phase
\begin{equation} \theta_x^2+\theta_y^2=\sigma^2, \end{equation}
where 
$\sigma ^2$ is solved with the dispersion relation, taking into account the variable depth and 
$\omega$. The rays, orthogonal to wavefronts, are searched by applying the method of characteristics to (4.7). We refer the reader to Johnson (Reference Johnson1997) and Zauderer (Reference Zauderer2006) for details of this method for solving the eikonal equation.
Figure 16 illustrates the comparison between the level curves of numerical solutions and the theoretically calculated rays. Clearly, for the wave profiles (
$\eta _1$ and 
$\eta _2$), the deformation of wavefronts in both modes corresponds well to the ray directions. Furthermore, wave profiles at different instants are presented to find the location of the maximal amplitude. The focusing point of wave energy is located in the vicinity of the position where the rays meet. It is found that the maximum wave elevation on the lower interface is situated before the one on the upper interface is located. A reasonable explanation is that internal waves on the lower interface interact directly with the bottom topography, and the boundary effect of the mound influences the position of the focusing point. Moreover, the topography effect is transferred to the upper layer by coupling between interfacial waves, which have a delay to some extent. The focusing point predicted by theoretical solutions for the mode 1 internal wave is 
$(7.72,5.5)$, and the average position where the wave profiles on both interfaces attain the maximal elevation at 
$t\approx 10$ is 
$(7.73,5.5)$. The relative error is apparently tiny. In these examples, with a non-trivial amplitude, our three-dimensional simulations capture the predicted dynamic features of interfacial waves well.
A simulation for a mode 2 internal wave is also shown in figure 17. In contrast to mode 1, the mode 2 wave on two interfaces exhibits opposite phases. Additionally, in figures 17(b) and 17(d) the envelopes of wave profiles at different instants also show opposite variations. Notably, the position where the amplitude of 
$\eta _1$ reaches its maximum is a local minimum of 
$\eta _2$. Considering the focusing of 
$\eta _1$, which is more direct, the focusing point is at 
$(5.5,8.3)$ while the corresponding local extremum of 
$\eta _2$ is at 
$(5.5,8.4)$. The theoretical focusing point position is 
$(5.5,8.33)$, close to our results. The delay phenomenon still exists but is less obvious.
Figure 17. A comparison between numerical solutions and theoretically predicted rays for a mode-2 incident internal wave. In the upper figures the rays are presented with dashed lines, and the projections of interfaces are plotted for 
$\eta _1$ (a) and 
$\eta _2$ (b). Wave profiles along the streamwise direction at 
$y=y_c=5.5$ are shown in the bottom figures for 
$t=13$ (solid lines), 
$t=14$ (dashed lines), 
$t=15$ (dotted lines) and 
$t=16$ (dash-dotted lines). The amplitude of 
$\eta _1$ reaches its maximum at 
$x\approx 8.3$.
4.2.2. Bathymetry in the strait
In this section we investigate a two-fluid system bounded above by a free surface, a limiting case of a three-layer fluid configuration, assuming the density of the top layer is negligible. The theoretical derivation and numerical procedure for this situation, which are akin to the three-layer case, are presented in Appendix A. Figure 18(a) displays the bathymetry of the Strait of Gibraltar within an 
$80 \times 80$ km horizontal area, sourced from the Smith and Sandwell 2-minute database. The data were downloaded from Smith and Sandwell 2-minute Database (2004), and the interested readers are referred to Marks & Smith (Reference Marks and Smith2006) for more global bathymetry data sets. This geographical location is a focal point of our investigation due to its unique bathymetric characteristics. The original latitude and longitude data were converted into a Cartesian coordinate system and interpolated with the bilinear method on the numerical grid. The depth of the upper-layer fluid is taken as 100 m. The area covered by land over water is treated with shallow-water conditions, where we take 100 m as the depth of the lower-layer fluid. Although this treatment is enough for accurate wave simulation in the area of the strait, some unphysical high-frequency components in the spanwise 
$y$ direction (orthogonal to the direction of wave propagation) gradually appear in the time evolution and need to be absorbed during calculation. We also remark that the largest depth in the computational domain is about 1200 m. Here 
$H_1=800$ m and 
$l=8000$ m are used for non-dimensionalization as reference characteristic lengths. The ratio of density 
$R=0.9$ is selected to be close to 1 to describe a small density jump.
Figure 18. Illustration of the real bathymetry and numerical simulation of the interfacial and free-surface waves passing through the Strait of Gibraltar: (a) bathymetry; (b) internal wave contours at 
$t=6$; (c) satellite image of internal waves in the Strait of Gibraltar; (d–f) successive snapshots for wave profiles at 
$t=1.2$, 
$t=3$ and 
$t=5.5$.
The initial perturbation is introduced as a Gaussian velocity potential along the interface in the streamwise 
$x$ direction centred around 
$x=2$, where the bottom topography height is approximately 0 (which means the lower-layer depth is close to 100 m). The initial velocity potentials and interface displacements are
\begin{equation} \phi_1(x,y,1)={\rm e}^{{-}36(x-x_c)^2}\quad\text{and} \quad\phi_2=\eta_1=\eta_2=0. \end{equation} This perturbation is initially uniform in the spanwise 
$y$ direction and then restricted to the strait using a filter function. The superposition of internal wave contours and varying bathymetry, calculated after six-time units with bathymetry, is presented in figure 18(b), alongside a satellite image of internal waves (see figure 18c). The wavefront exhibits a bow-shaped pattern upon exiting the strait, with diffraction occurring in both the southern and northern orientations. Local peaks in wave elevation form a distribution curve originating from the strait's centre and extending towards the southeast. From figure 18(b), the deep and shallow areas can be clearly distinguished. In the shallow area, the wavenumber increases and the phase velocity decreases. In contrast, the wavefronts propagate faster in the deep area, which leads to obvious inflection along the deepest region. Applying two lines of references, these observations are well aligned with numerical results in figures 18(d–f) and 19(a) and the satellite image.
Figure 19. Free-surface projections at 
$t=6$ calculated with (a) bathymetry in the Strait of Gibraltar and (b) a flat bottom. (c–e) Cross-sections at different positions: interfacial waves with bathymetry (solid line), free-surface waves with bathymetry (dashed line), interfacial waves with a flat bottom (square mark) and free-surface waves with a flat bottom (dashed line with circular mark). (c) Wave profiles at 
$y=5$ and the topography profile at the same position marked with a red dotted line; (d) wave profiles at 
$y=6$; (e) wave profiles with bathymetry at 
$x=6.8$ and with a flat bottom at 
$x=6.3$.
Figure 19 allows us to analyse the topographic impact by comparing it to wave fields in the flat bottom case. The Strait of Gibraltar's bathymetry restricts the width of wave development and leads to asymmetry in wave propagation. Energy becomes more concentrated in a narrow space near wave peaks. This effect is further explored in figure 19(c–e), where wave profiles at different positions are extracted and superimposed for comparison. Mode 2 internal and free-surface waves (i.e. the free surface and interface are of the opposite phases) exist in the vicinity of the initial perturbance. The region dominated by mode 2 waves grows in time, and the amplitude of a free-surface wave is much smaller than that of the associated internal wave. Further from the centre of the initial Gaussian disturbance in the velocity potential on the interface, waves gradually transform into mode 1, and the free-surface wave admits an amplitude in the same order as the interfacial wave. A flat platform that is prolonged with the increase of time forms between two regions. These phenomena are also recovered in the flat bottom case. However, the rising bottom topography in the second region leads to a smaller wavelength (see figure 19c). With the observing point moving from the spanwise centre to the north in figure 19(d), the amplitude decrease is more rapid in the strait case. The asymmetric character in the spanwise direction is presented in figure 19(e).
5. Conclusions
We have developed a numerical method featuring DtN operators, enabling investigations of multimode internal waves passing over significant bottom topography. This approach operates within the framework of linear approximation, eliminating limitations on the height or smoothness of bottom topography. Moreover, our numerical model is adaptable to three-dimensional multi-layer systems. We employ the Galerkin approximation to save computational costs while accurately capturing typical three-dimensional effects.
We initially validated our numerical method by applying it to two-dimensional linear wave shoaling cases. In a three-fluid system we explored the influence of various parameters on wave behaviour, including depth and density ratios. Multiple simulations with different parameter settings revealed changing trends in wave amplitudes. Notably, we identified a critical solution that yielded intriguing results. Additionally, within a three-layer system, we analysed the excitation of waves of different modes induced by bottom topography. Our numerical simulations of wavepackets passing over locally confined bottom obstacles demonstrated that the obstacle's slope plays a pivotal role in exciting mode 2 internal waves, aligning with the findings of Liu et al. (Reference Liu, Grimshaw and Johnson2019). We also examined the effects of rapidly varying bottom topography. By considering internal wavepackets with a fixed typical wavelength and changing the fundamental wavelength of bottom topography, we observed classic and novel phenomena, including Bragg resonance, mode 2 internal wave excitation, special harmonics resonance and the homogenization effect of quickly varying bottom topography.
For three-dimensional simulations, our numerical procedure, coupled with the Galerkin method, accurately replicated the deformation of wavefronts and the position of a focusing point in the presence of a Luneberg lens mound compared with theoretical values. We also identified a subtle delay in how upper-layer internal waves respond to topography effects. We used real topographic data from the Strait of Gibraltar for simulation in a two-and-a-half-layer system featuring an interface and a free surface. Our results successfully replicated the propagation directions of waves compared with satellite imageries. Furthermore, our study provided insights into the coupled evolution between internal and free-surface waves, expanding upon previous research.
Finally, we remark that to include both wave nonlinearity and significant bottom topography in the three-layer fluid system to understand the dynamics of moderate- or large-amplitude multimode internal waves (for example, generation of high-mode solitary-like internal waves due to the resonance effect of bottom topography, transitions between nonlinear internal waves of different modes through bottom topography, fission of high-mode solitary-like internal waves propagating over sharply varying topography, etc.) merits a thorough investigation, which we leave for future work.
Funding
This work was supported by the National Science Foundation for Distinguished Young Scholars (no. 12325207) and the Key Program of the National Natural Science Foundation of China (no. 12132018).
Declaration of interest
The authors report no conflict of interest.
Appendix A. Two-layer configuration in the presence of a free surface
A.1. Governing equations and linear dispersion relation
A sketch of a two-fluid system bounded by a free surface on top (sometimes termed the two-and-a-half-layer system) is given in figure 20, and the corresponding governing equations are
\begin{equation} \left.\begin{array}{c} \phi_{1xx}+\phi_{1yy}+\phi_{1zz}=0 \quad \text{for } b(x,y)< z< H_1+\eta_1(x,y,t),\\ \phi_{2xx}+\phi_{2yy}+\phi_{2zz}=0 \quad \text{for } H_1+ \eta_1(x,y,t)< z< H_1+H_2+\eta_2(x,y,t). \end{array} \right\} \end{equation}The impermeability boundary condition at the bottom writes
\begin{equation} \phi_{1z}-\phi_{1x}b_x-\phi_{1y}b_y=0 \quad \text{on } z=b(x,y). \end{equation}The kinematic and dynamic conditions on the interface and free surface read
\begin{equation} \left.\begin{array}{c} \eta_{1t}+\eta_{1x}\phi_{1x}+\eta_{1y}\phi_{1y}-\phi_{1z}=0 \quad \text{on } z=H_1+\eta_1(x,y,t), \\ \eta_{1t}+\eta_{1x}\phi_{2x}+\eta_{1y}\phi_{2y}-\phi_{2z}=0 \quad \text{on } z=H_1+\eta_1(x,y,t), \\ \eta_{2t}+\eta_{2x}\phi_{2x}+\eta_{2y}\phi_{2y}-\phi_{2z}=0 \quad \text{on } z=H_1+H_2+\eta_2(x,y,t), \end{array} \right\} \end{equation}
\begin{equation} \left.\begin{array}{c} \begin{aligned} & \phi_{1t}-R\phi_{2t}+g(1-R)\eta_1+\tfrac{1}{2}(\phi_{1x}^2 +\phi_{1y}^2+\phi_{1z}^2-R\phi_{2x}^2-R\phi_{2y}^2-R\phi_{2z}^2)=0 \\ & \quad \text{on } z=H_1+\eta_1(x,y,t), \end{aligned}\\ \phi_{2t}+g\eta_2+\tfrac{1}{2}(\phi_{2x}^2+\phi_{2y}^2+\phi_{2z}^2)=0 \quad \text{on } z=H_1+H_2+\eta_2(x,y,t). \end{array} \right\} \end{equation}
Figure 20. Sketch of a two-layer density stratified system separated by an interface and bounded above by a free surface.
Applying the scaling (2.5a–g) to linearize the whole system, the equations to be numerically solved write
\begin{equation} \left.\begin{array}{c} \mu^2(\phi_{1xx}+\phi_{1yy})+\phi_{1zz}=0 \quad \text{for } b< z<1,\\ \mu^2(\phi_{2xx}+\phi_{2yy})+\phi_{2zz}=0 \quad \text{for } 1< z<1+\varTheta,\\ \phi_{1z}-\mu^2(\phi_{1x}b_x+\phi_{1y}b_y)=0 \quad \text{on } z=b,\\ \eta_{1t}-\dfrac{1}{\mu^2}\phi_{1z}=\eta_{1t}- \dfrac{1}{\mu^2}\phi_{2z}=0 \quad \text{on } z=1,\\ \eta_{2t}-\dfrac{1}{\mu^2}\phi_{2z}=0 \quad \text{on } z=1+\varTheta,\\ \phi_{1t}-R\phi_{2t}+(1-R)\eta_1=0 \quad \text{on } z=1,\\ \phi_{2t}+\eta_2=0 \quad \text{on } z=1+\varTheta, \end{array} \right\} \end{equation}
where 
$\varTheta =H_2/H_1$ and 
$R<1$ is the upper-to-lower-layer density ratio. Then, performing a similar derivation introduced in § 2.2, one obtains
\begin{align} &\left(\frac{\mu\omega_0^2}{k^2}\right)^2\left(\frac{k^2}{R}+k^2\tanh(\mu k)\tanh(\mu k\varTheta)\right)-\frac{\mu\omega_0^2}{k^2}\left(\frac{k}{R}\tanh(\mu k)+\frac{k}{R}\tanh(\mu k\varTheta)\right)\nonumber\\ &\quad +\left(\frac{1}{R}-1\right)\tanh(\mu k)\tanh(\mu k\varTheta)=0, \end{align}
which represents the dispersion relation of a two-and-a-half-layer system. It is worth noting that the coefficients in this quadratic equation are the same as those obtained directly from (2.19) by taking the limits 
$R_2\to 0$ and 
$\varTheta _2\to +\infty$. Moreover, with 
$(k, \mu, \varTheta )>0$ and 
$0< R<1$, both solutions to (A6) are positive, indicating that the whole space of system parameters remains feasible with real phase velocity. The ratio of wave displacement between the free surface and the interface can be expressed as
\begin{equation} \frac{a_2}{a_1}=\frac{\cosh(\mu k\varTheta)}{R}\left[\frac{\mu\omega_0^2}{k}\left(\frac{1}{\tanh(\mu k)}+R\tanh(\mu k\varTheta)\right)+R-1\right]. \end{equation}A.2. Numerical formulation
We first divide the linear system into two elliptic boundary value problems and one time-evolution system by introducing the DtN operators. Two boundary value problems are
\begin{equation} \left.\begin{array}{c} \mu^2(\phi_{1xx}+\phi_{1yy})+\phi_{1zz}=0 \quad \text{for } b< z<1, \\ \phi_{1z}-\mu^2(\phi_{1x}b_x+\phi_{1y}b_y)=0 \quad \text{on } z=b, \\ \phi_1=q_1\quad \text{on } z=1, \end{array} \right\} \end{equation}and
\begin{equation} \left.\begin{array}{c} \mu^2(\phi_{2xx}+\phi_{2yy})+\phi_{1zz}=0 \quad \text{for } 1< z<1+\varTheta, \\ \phi_2=q_{22} \quad \text{on } z=1+\varTheta, \\ \phi_2=q_{12} \quad \text{on } z=1, \end{array} \right\} \end{equation}
where we introduce 
$q_1(\boldsymbol {x},t)=\phi _1(\boldsymbol {x},1+\epsilon \eta _1(\boldsymbol {x},t),t)$, 
$q_{12}(\boldsymbol {x},t)=\phi _2(\boldsymbol {x},1+\epsilon \eta _1(\boldsymbol {x},t),t)$ and 
$q_{22}(\boldsymbol {x},t)=\phi _2(\boldsymbol {x},1+\varTheta +\epsilon \eta _2(\boldsymbol {x},t),t)$ as three new variables. We remark that with the kinematic boundary condition on the interface
\begin{equation} \phi_{1z}(\boldsymbol{x},z=1,t)=\phi_{2z}(\boldsymbol{x},z=1,t), \end{equation}
one can obtain 
$q_{12}=\mathcal {H}[q_1, q_{22}]$. Furthermore, for the boundary value problem in the upper layer, we notice that unlike in the lower layer, the problem has two Dirichlet boundary conditions rather than one Dirichlet and one Neumann condition. The time evolution is written as
\begin{equation} \left.\begin{array}{c} \eta_{1t}=\dfrac{1}{\mu^2}\mathcal{G}_1[q_1], \\ \eta_{2t}=\dfrac{1}{\mu^2}\mathcal{G}_2[q_{22},q_{1}], \\ q_{1t}-Rq_{12t}=(R-1)\eta_1,\\ q_{22t}={-}\eta_2, \end{array} \right\} \end{equation}
where 
$\mathcal {G}_1$ and 
$\mathcal {G}_2$ are the DtN operators defined as 
$\mathcal {G}_1[q_1]=\phi _{1z}(\boldsymbol {x},z=1,t)=\phi _{2z}(\boldsymbol {x},z=1,t)$ and 
$\mathcal {G}_2[q_{22},q_1]=\phi _{2z}(\boldsymbol {x},z=1+\varTheta,t)$, respectively.
Applying the Fourier transform to the elliptic problems, the expression of the velocity potential in the lower layer has the same form as (3.13), namely
\begin{equation} \phi_1(x,y,z)=\widehat{q_1}(\boldsymbol{0})+ \sum_{\boldsymbol{k}\in\varLambda^*}\, {\rm e}^{{\rm i}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\left[\widehat{q_1}(\boldsymbol{k})\frac{\cosh(\mu kz)}{\cosh(\mu k)}+X_{\boldsymbol{k}}\frac{\sinh(\mu k(z-1))}{k\cosh^2(\mu k)}\right], \end{equation}and the solution of topographic coefficients follows the same procedure described in § 3.2. The solution of the second boundary value problem is
\begin{equation} \left.\begin{array}{c} \displaystyle \phi_2(x,y,z)=\sum_{\boldsymbol{k}\in\varLambda^*}\, {\rm e}^{{\rm i}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\left[\widehat{q_{22}}(\boldsymbol{k})\dfrac{\sinh(\mu k(z-1))}{\sinh(\mu k\varTheta)}+\widehat{q_{21}}(\boldsymbol{k})\dfrac{\sinh(\mu k(1+\varTheta-z))}{\sinh(\mu k \varTheta)}\right],\\ \displaystyle \phi_{2z}(x,y,z=1)=\sum_{\boldsymbol{k}\in\varLambda^*}\, {\rm e}^{{\rm i}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\left[\widehat{q_{22}}(\boldsymbol{k})\dfrac{\mu k}{\sinh(\mu k \varTheta)}-\widehat{q_{12}}(\boldsymbol{k})\dfrac{\mu k}{\tanh(\mu k \varTheta)}\right],\\ \displaystyle \phi_{2z}(x,y,z=1+\varTheta)=\sum_{\boldsymbol{k}\in\varLambda^*}\, {\rm e}^{{\rm i}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}\left[\widehat{q_{22}}(\boldsymbol{k})\dfrac{\mu k}{\tanh(\mu k \varTheta)}-\widehat{q_{12}}(\boldsymbol{k})\dfrac{\mu k}{\sinh(\mu k \varTheta)}\right], \end{array}\right\} \end{equation}
where we suppose, without loss of generality, 
$\widehat {q_{22}}(\boldsymbol {0})=\widehat {q_{12}}(\boldsymbol {0}) =\widehat {\phi _2}(\boldsymbol {0})=0$. With the orthogonality of the Fourier primary functions, the matching condition on the interface requires
\begin{equation} \widehat{q_{12}}(\boldsymbol{k})={-}\tanh(\mu k)\tanh(\mu k\varTheta)\widehat{q_1}(\boldsymbol{k})-\frac{\tanh(\mu k \varTheta)}{k\cosh^2(\mu k)}X_{\boldsymbol{k}}+\frac{1}{\cosh(\mu k\varTheta)}\widehat{q_{22}}(\boldsymbol{k}), \end{equation}which signifies that
\begin{equation} \widehat{q_{12}}_t(\boldsymbol{k})={-}\tanh(\mu k) \tanh(\mu k\varTheta)\widehat{q_1}_t(\boldsymbol{k})- \frac{\tanh(\mu k\varTheta)}{k\cosh^2(\mu k)}X_{\boldsymbol{k}t}+ \frac{1}{\cosh(\mu k \varTheta)}\widehat{q_{22}}_t(\boldsymbol{k}). \end{equation}Finally, the time-evolution problem writes
\begin{gather} \left.\begin{array}{c} \displaystyle \widehat{\eta_1}_t(\boldsymbol{k})=\dfrac{1}{\mu^2} \widehat{\mathcal{G}_1[q_1]}(\boldsymbol{k})= \dfrac{1}{\mu^2}\left[\widehat{q_1}(\boldsymbol{k})\mu k\tanh(\mu k)+ X_{\boldsymbol{k}}\mu\dfrac{1}{\cosh^2(\mu k)}\right],\\ \begin{aligned} & \widehat{\eta_2}_t(\boldsymbol{k})=\dfrac{1}{\mu^2} \widehat{\mathcal{G}_2[q_{22},q_1]}(\boldsymbol{k})= \dfrac{1}{\mu^2}\widehat{q_{22}}(\boldsymbol{k})\left(\dfrac{\mu k}{\tanh(\mu k \varTheta)}-\dfrac{\mu k }{\cosh(\mu k \varTheta)\sinh(\mu k \varTheta)}\right)\\ & \quad +\widehat{q_1}(\boldsymbol{k})\dfrac{\mu k\tanh(\mu k)}{\cosh(\mu k \varTheta)}+X_{\boldsymbol{k}}\dfrac{\mu}{\cosh^2(\mu k)\cosh(\mu k \varTheta)}, \end{aligned} \end{array} \right\} \end{gather}
\begin{gather} \left.\begin{aligned} &\widehat{q_1}_t(\boldsymbol{k})(1+R\tanh(\mu k)\tanh(\mu k\varTheta))+\dfrac{R\tanh(\mu k\varTheta)}{k\cosh^2(\mu k)}X_{\boldsymbol{k}t}\\ &\quad =(R-1)\widehat{\eta_1}(\boldsymbol{k})-\dfrac{R\tanh(\mu k\varTheta)}{\sinh(\mu k\varTheta)}\widehat{\eta_2},\\ &\widehat{q_2}_t(\boldsymbol{k})={-}\widehat{\eta_2}. \end{aligned}\right\} \end{gather}
