2.1. Basic equations

Consider a two-layer system composed of three-dimensional incompressible, inviscid and homogeneous fluids, with a lighter fluid on top of a heavier one (see the schematic in figure 2). Two fluids are assumed to be immiscible and separated by a sharp interface $z=\eta (x,y,t)$ , where $x$ and $y$ are horizontal coordinates, and the $z$ -axis points upwards with $z=0$ the undisturbed interface. The upper layer is semi-infinite, and the lower layer is bounded below by a bottom topography $z=-h+b(x,y,t)$ varying in space and time, where $b(x,y,t)$ is a given function. It can be used to model a non-uniform-in-time background current over a bottom obstacle, or underwater landslides interacting with density-stratified fluids. The fluid density in each layer is assumed to be constant, denoted by $\rho ^+$ and $\rho ^-$ for the upper and lower layers, respectively, with $\rho ^+\lt \rho ^-$ representing a stable density configuration. Flows in their respective domains are assumed to be irrotational, therefore the velocity potential in each layer satisfies the Laplace equation, namely,

(2.1) \begin{align} &\varDelta \phi ^-+\phi _{zz}^-=0\quad \text{for } \quad {-h}+b(x,y,t)\lt z\lt \eta (x,y,t),\nonumber \\ &\varDelta \phi ^++\phi _{zz}^+=0\quad \text{for } \quad z\gt \eta (x,y,t), \end{align}

where $\varDelta =\partial _{xx}+\partial _{yy}$ is the two-dimensional horizontal Laplace operator. On the interface $z=\eta (x,y,t)$ , the nonlinear kinematic and dynamic boundary conditions read

(2.2) \begin{equation} \eta _t={\phi _z^-} -\boldsymbol{\nabla }{\phi ^-}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta ={\phi ^+_z}-\boldsymbol{\nabla }{\phi ^+}\boldsymbol{\cdot }\boldsymbol{\nabla }\eta , \end{equation}

(2.3) \begin{equation} \rho ^-\left [\phi ^-_t+\frac 12|\boldsymbol{\nabla }\phi ^-|^2+\frac 12(\phi _z^-)^2+g\eta \right ]=\rho ^+\left [\phi ^+_t+\frac 12|\boldsymbol{\nabla }\phi ^+|^2+\frac 12(\phi _z^+)^2+g\eta \right ]\!, \end{equation}

where $\boldsymbol{\nabla}$ and $\boldsymbol{\nabla }\boldsymbol{\cdot}$ are the horizontal gradient and divergence operators, respectively, and $g$ is the acceleration due to gravity. The kinematic boundary condition (2.2) implies that the normal velocity is continuous across the interface, namely,

(2.4) \begin{eqnarray} \big (\phi ^-_x,\phi _y^-,\phi _z^-\big )^{\rm T} \boldsymbol{\cdot }\boldsymbol {n}=\big (\phi ^+_x,\phi ^+_y,\phi ^+_z\big )^{\rm T} \boldsymbol{\cdot }\boldsymbol {n}, \end{eqnarray}

where $\boldsymbol{n}=(-\eta _x,-\eta _y,1)^{\rm T} \big /\sqrt {1+|\boldsymbol{\nabla }\eta |^2}$ is the unit normal vector pointing upwards from the interface. In addition, the boundary conditions

(2.5) \begin{eqnarray} \phi ^+_z\rightarrow 0\quad \text{at } \quad z\rightarrow \infty, \end{eqnarray}

(2.6) \begin{eqnarray} \phi ^-_z-\boldsymbol{\nabla }b\boldsymbol{\cdot }\boldsymbol{\nabla }\phi ^-=b_t\quad \text{at } \quad z=-h+b(x,y,t), \end{eqnarray}

close the system.

Figure 2.

Sketch of the problem and coordinate system.

2.2. Ablowitz–Fokas–Musslimani formulation

Ablowitz, Fokas & Musslimani (Reference Ablowitz, Fokas and Musslimani2006) proposed an explicit non-local formulation for the classical water wave problem in both two and three dimensions. It was later generalised to include the bottom topography by Fokas & Nachbin (Reference Haut and Ablowitz2012), Curtis & Shen (Reference Djordjevic and Redekopp2015) and Andrade & Nachbin (Reference Andrade and Nachbin2018), and to study interfacial waves by Haut & Ablowitz (Reference Havelock2009). For completeness, we briefly introduce the Ablowitz–Fokas–Musslimani formulation. It is straightforward to verify that the identity

(2.7) \begin{equation} \left (\phi _z^\pm \psi _x+\phi _x^\pm \psi _z\right )_x+\big (\phi _z^\pm \psi _y+\phi _y^\pm \psi _z\big)_y+\big(\phi _z^\pm \psi _z-\phi _x^\pm \psi _x-\phi _y^\pm \psi _y\big)_z=0 \end{equation}

holds for an arbitrary harmonic function $\psi$ . Integrating (2.7) over the lower layer and applying the divergence theorem, one obtains

(2.8) \begin{align} \begin{aligned} 0=&\iint \big [b_x\big (\phi _z^-\psi _x+\phi _x^-\psi _z\big )+b_y\big (\phi ^-_z\psi _y+\phi ^-_y\psi _z\big)-\big (\phi ^-_z\psi _z-\boldsymbol{\nabla }\phi ^-\boldsymbol{\cdot }\boldsymbol{\nabla }\psi \big)\big ]_{z=-h+b} \mathrm{d}\boldsymbol {r}\\ &+\iint \big [-\eta _x\big (\phi _z^-\psi _x+\phi _x^-\psi _z\big )-\eta _y\big (\phi ^-_z\psi _y+\phi ^-_y\psi _z\big)+\big (\phi ^-_z\psi _z-\boldsymbol{\nabla }\phi ^-\boldsymbol{\cdot }\boldsymbol{\nabla }\psi \big )\big ]_{z=\eta } \mathrm{d}\boldsymbol {r},\\[3pt] \end{aligned} \end{align}

where $\boldsymbol {r}=(x,y)^{\rm T}$ is a vector in the horizontal plane. Substituting $\psi =\textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}+kz}$ into (2.8) yields

(2.9) \begin{align} 0=\iint \textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}+k(\eta +h)}\big [k\eta _t-\text{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol{\nabla }\varPhi ^-\big ]\,\mathrm{d}\boldsymbol {r}+\iint \textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}+kb}\big[-kb_t+\text{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol{\nabla }\varPhi ^b\big ]\,\mathrm{d}\boldsymbol {r}, \end{align}

where $\varPhi ^-=\phi ^-(x,y,\eta ,t)$ and $\varPhi ^b=\phi ^-(x,y,b,t)$ are interfacial and bottom velocity potentials, respectively, $\boldsymbol {k}=(k_x,k_y)$ is the wavenumber vector, and $k=|\boldsymbol {k}|=\sqrt {k_x^2+k_y^2}$ is its modulus. Since $\textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}-kz}$ is also a solution to the Laplace equation, one then obtains

(2.10) \begin{align} 0=\iint \textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}-k(\eta +h)}\big [-k\eta _t-\text{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol{\nabla }\varPhi ^-\big ]\, \mathrm{d}\boldsymbol {r}+\iint \textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}-kb}\big[kb_t+\text{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol{\nabla }\varPhi ^b\big ]\, \mathrm{d}\boldsymbol {r}. \end{align}

Adding and subtracting (2.9) and (2.10) yields

(2.11) \begin{equation} \begin{aligned} 0={}&\iint \textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}}\left [k\sinh (k(\eta +h))\,\eta _t-{\rm i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol{\nabla }\varPhi ^-\cosh (k(\eta +h))\right ]\, \mathrm{d}\boldsymbol {r}\\ &{}+\iint \textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}}\big [-k\sinh (kb)\,b_t+{\rm i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol{\nabla }\varPhi ^b\cosh (kb)\big ]\,\mathrm{d}\boldsymbol {r} \end{aligned} \end{equation}

and

(2.12) \begin{eqnarray} \begin{aligned} 0={}&\iint \textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}}\left [\cosh (k(\eta +h))\,\eta _t- \text{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol{\nabla }\varPhi ^{-\frac {\sinh (k(\eta +h))}{k}}\right ]\mathrm{d}\boldsymbol {r}\\[3pt] &{}+\iint \textrm{e}^{\text{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}}\left [-\cosh (kb)\,b_t+\text{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol{\nabla }\varPhi ^b\frac {\sinh (kb)}{k}\right ]\mathrm{d}\boldsymbol {r}. \end{aligned} \end{eqnarray}

In the same vein, one can obtain the global relation for the upper layer

(2.13) \begin{equation} 0=\iint \textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}}\left [k\eta _t+\text{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol{\nabla }\varPhi ^+\right ]\textrm{e}^{-k\eta }\,\mathrm{d}\boldsymbol {r}, \end{equation}

where $\varPhi ^+=\phi ^+(x,y,\eta ,t)$ . Finally, if we denote the density ratio by $R=\rho ^+/\rho ^-\lt 1$ , then the dynamic boundary condition on the interface can be rewritten as

(2.14) \begin{equation} \begin{aligned} 0={}&\left ({\varPhi ^-} -{R\varPhi ^+}\right )_t+(1-R)g\eta +\frac 12\left |\boldsymbol{\nabla }{\varPhi ^-}\right |^2-\frac {R}{2}\left |\boldsymbol{\nabla }{\varPhi ^+}\right |^2\\[5pt] &{}-\frac 12\frac {\big(\eta _t+\boldsymbol{\nabla }\eta \boldsymbol{\cdot }\boldsymbol{\nabla }{\varPhi ^-}\big)^2}{1+|\boldsymbol{\nabla }\eta |^2}+\frac {R}{2}\frac {\big(\eta _t+\boldsymbol{\nabla }\eta \boldsymbol{\cdot }\boldsymbol{\nabla }{\varPhi ^+}\big)^2}{1+|\boldsymbol{\nabla }\eta |^2}. \end{aligned} \end{equation}

In the remainder of the paper, we consider the situation where the typical wavelength is large compared to the mean depth of the lower-layer fluid. Under such circumstances, the linear dispersion relation and its approximation can be expressed as

(2.15) \begin{equation} \omega ^2=\frac {({\rho ^-}-{\rho ^+})gk\tanh (kh)}{{\rho ^-}+{\rho ^+}\tanh (kh)}\approx (1-R)ghk^2 \end{equation}

for $kh\ll 1$ , which gives a reference speed for long-wave models, namely $\sqrt {(1-R)gh}$ . In the subsequent analyses, two types of model equations will be derived: fully nonlinear models (i.e. no assumption is made about the wave amplitude), and weakly nonlinear models (i.e. with a suitable smallness assumption on the wave amplitude).

3.1. Vertically averaged variables

The system (3.8), (3.13) can also be derived from vertically averaged quantities without any irrotationality assumption for the lower layer. Let $(\boldsymbol {u}^-,w^-)=(u^-,v^-,w^-)$ be the velocity components of the lower-layer fluid, and let $p^-$ be the pressure. These unknown variables satisfy the incompressible and inviscid Navier–Stokes equations

(3.16) \begin{equation} \left \{ \begin{aligned} &\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol {u}^-}+{\partial _zw^-}=0,\\ &{\boldsymbol {u}^-_t}+({\boldsymbol {u}^-}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol {u}^-}+{w^-}\partial _z{\boldsymbol {u}^-}={-\boldsymbol{\nabla }p^-},\\ &\mu ^2\left ({w^-_t}+{\boldsymbol {u}^-}\boldsymbol{\cdot }{\boldsymbol{\nabla }w^-}+{w^-}{\partial _z}{w^-}\right )=-{\partial _z}{p^-}-\dfrac{1}{1-R}. \end{aligned}\right . \end{equation}

The above Navier–Stokes equations have been non-dimensionalised by assuming

(3.17) \begin{align} \begin{aligned} &x,y\sim \lambda ,\quad z\sim h,\quad t\sim \frac {\lambda }{\sqrt {(1-R)gh}},\\ &\boldsymbol {u}^-\sim \sqrt {(1-R)gh},\quad w^-\sim \frac {h}{\lambda }\sqrt {(1-R)gh},\quad p^-\sim (1-R)\rho ^-gh. \end{aligned} \end{align}

The last equation of (3.16) implies that

(3.18) \begin{equation} p^-(x,y,z,t)=-\frac {z}{1-R}+P(x,y,t)+O(\mu ^2). \end{equation}

Integrating the equation of conservation of mass in the vertical direction from $z=-1+b(x,y,t)$ to $z=\eta (x,y,t)$ gives

(3.19) \begin{equation} (\eta -b)_t+\boldsymbol{\nabla }\boldsymbol{\cdot }\left [(1+\eta -b)\overline {\boldsymbol {u}^-}\right ]=0, \end{equation}

where the kinematic boundary conditions

(3.20) \begin{align} &w^-=\eta _t+\boldsymbol {u}^-\boldsymbol{\cdot }\boldsymbol{\nabla }\eta\quad \text{at } \quad z=\eta (x,y,t),\nonumber \\ &w^-=b_t+\boldsymbol {u}^-\boldsymbol{\cdot }\boldsymbol{\nabla }b\quad \text{at } \quad z=-1+b(x,y,t), \end{align}

are used to replace $w^-$ whenever necessary, and

(3.21) \begin{align} \overline {f}(x,y,t)=\frac {1}{1+\eta -b}\int _{-1+b}^\eta f(x,y,z,t)\,\mathrm{d}z. \end{align}

Integrating the horizontal momentum equations yields

(3.22) \begin{equation} \left [(1+\eta -b)\overline {\boldsymbol {u}^-}\right ]_t+\boldsymbol{\nabla }\boldsymbol{\cdot }\left [(1+\eta -b)\overline {\boldsymbol {u}^-}\;\overline {\boldsymbol {u}^-}\right ]=-(1+\eta -b)\,\boldsymbol{\nabla }\overline {P}+O(\mu ^2). \end{equation}

To close the system, the pressure $P$ needs to be determined. If the flow in the upper layer is assumed to be irrotational, then the Bernoulli equation in the upper layer gives the expression of the pressure at the interface:

(3.23) \begin{equation} p^+=-R\left [\varPhi ^+_t+\frac 12|\boldsymbol{\nabla }\varPhi ^+|^2+\frac {\eta }{1-R}-\frac {\mu ^2}{2}\frac {(\eta _t+\boldsymbol{\nabla }\eta \boldsymbol{\cdot }\boldsymbol{\nabla }\varPhi ^+)^2}{1+\mu ^2\,|\boldsymbol{\nabla }\eta |^2}\right ]\!. \end{equation}

Due to (3.9) and the continuity of pressure across the interface, one obtains

(3.24) \begin{equation} P=p^-+\frac {\eta }{1-R}=p^++\frac {\eta }{1-R}=\eta +\mu \textit{RG}_0^{-1}\eta _{tt}+O(\mu ^2). \end{equation}

Substituting (3.24) into (3.22), and neglecting terms of order $\mu ^2$ , one obtains, after some algebra,

(3.25) \begin{equation} \overline {\boldsymbol {u}^-}_t+\boldsymbol{\nabla }\eta +\overline {\boldsymbol {u}^-}\boldsymbol{\cdot }\boldsymbol{\nabla }\overline {\boldsymbol {u}^-}+\mu R\,\boldsymbol{\nabla }G_0^{-1}\eta _{tt}=0. \end{equation}

It is evident that the system (3.19), (3.25) is precisely the same as (3.8), (3.13).

Finally, it is worth mentioning that when both fluid layers are shallow compared to the typical wavelength, the Miyata–Choi–Camassa system can be further regularised by using the velocities at the top and bottom rigid walls as unknowns, as opposed to their vertical averages (for details, see Choi, Barros & Jo Reference Choi and Camassa2009; Barros & Choi Reference Barros and Choi2013). However, this method seems not applicable for the ‘shallow–deep’ configuration upon noticing $\varPhi ^b-\varPhi ^-=O (\mu ^2 )$ , so the modification cannot enter the system.

3.2. Radially symmetric model

The surface signature of ocean internal waves typically features curvilinear fronts, a ubiquitous phenomenon visible in synthetic aperture radar imagery from satellites. These wave fronts look like annular waves, and are usually approximated as solutions to the cylindrical KdV equation (see e.g. Khusnutdinova & Zhang Reference Lindzen2016). The curvilinear internal wave fronts also occur in the low atmosphere. A typical example of arched cloud formations that indicate atmospheric internal gravity waves is shown in figure 3 (Yiğit & Medvedev Reference Yiğit and Medvedev2019). Motivated by these recent observations, we write (3.8), (3.11) in polar coordinates. Introducing the coordinate transformation $x=r\cos \theta$ and $y=r\sin \theta$ , one obtains

(3.26) \begin{equation} \frac {\partial }{\partial x}=\cos \theta\, \frac {\partial }{\partial r}-\frac {\sin \theta }{r}\frac {\partial }{\partial \theta },\quad \frac {\partial }{\partial y}=\sin \theta\, \frac {\partial }{\partial r}+\frac {\cos \theta }{r}\frac {\partial }{\partial \theta }, \end{equation}

which implies

(3.27) \begin{equation} \varDelta =\frac {\partial ^2}{\partial r^2}+\frac {1}{r}\frac {\partial }{\partial r}+\frac {1}{r^2}\frac {\partial ^2}{\partial \theta ^2},\quad \boldsymbol{\nabla }f_1\boldsymbol{\cdot }\boldsymbol{\nabla }f_2=\frac {\partial f_1}{\partial r}\frac {\partial f_2}{\partial r}+\frac {\partial f_1}{\partial \theta }\frac {\partial f_2}{\partial \theta }. \end{equation}

For axisymmetric motions (i.e. $\partial _\theta =0$ ), we substitute (3.27) into (3.8), (3.11) to obtain the cylindrical MCC-BO equations

(3.28) \begin{equation} \left \{ \begin{aligned} &\eta _t-b_t+(\eta -b)_r\varPhi ^-_r+(1+\eta -b)\left (\varPhi ^-_{rr}+\varPhi ^-_r/r\right )=0,\\ &\varPhi ^-_t+\eta +\frac 12\left (\varPhi ^-_r\right )^2+\mu \textit{RG}_0^{-1}\eta _{tt}=0. \end{aligned}\right . \end{equation}

The pseudo-differential operator $G_0^{-1}$ in the axisymmetric case is defined in the following way. Since by definition

(3.29) \begin{equation} \mathcal{F}\left [G_0^{-1}\eta _{tt}\right ]=\frac {1}{k}\mathcal{F}\left [{\eta }_{tt}\right ]\Longrightarrow\quad G_0^{-1}\eta _{tt}=\mathcal{F}^{-1}\left [\frac {1}{k}\right ]*\eta _{tt}\,, \end{equation}

where the asterisk denotes convolution, $\mathcal{F}$ stands for the Fourier transform (which is termed the Hankel transform of order zero in the polar coordinates), and $\mathcal{F}^{-1}$ represents its inverse, it is not difficult to verify that

(3.30) \begin{equation} \mathcal{F}^{-1}\left [\frac {1}{k}\right ]=\left (\frac {1}{2\pi }\right )^2\iint \frac {\textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}}}{k}\,\mathrm{d}\boldsymbol {k}=\frac {1}{2\pi }\int _0^\infty J_0(kr)\,\mathrm{d}k=\frac {1}{2\pi r}, \end{equation}

where $r=|\boldsymbol {r}|$ , and $J_0$ is the zeroth-order Bessel function of the first kind. Therefore, for the special case of circularly symmetric functions that have no angular dependence, $\eta (\boldsymbol {r},t)$ , we have

(3.31) \begin{equation} \begin{aligned} G_0^{-1}\eta _{tt}=&\iint \eta _{tt}(\boldsymbol {r}',t)\;\left (\frac {1}{2\pi }\right )^2\iint \frac {\textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }(\boldsymbol {r}-\boldsymbol {r}')}}{k}\,\mathrm{d}\boldsymbol {k}\,\mathrm{d}\boldsymbol {r}'\\[3pt] =&\int _0^\infty r'\,\eta _{tt}(r',t)\int _0^\infty J_0(kr)\,J_0(kr')\,\mathrm{d}k\mathrm{d}r', \end{aligned} \end{equation}

where $\boldsymbol {r}'=(x',y')$ and $r'=|\boldsymbol {r}'|$ . It is noted that the inner integral can be simplified as

(3.32) \begin{equation} \int _0^\infty J_0(kr)\,J_0(kr')\,\mathrm{d}k=\begin{cases} \dfrac {2}{\pi r}\mathcal{K}\left (\dfrac {r'}{r}\right )&\text{for }\quad r\gt r',\\[11pt] \dfrac {2}{\pi r'}\mathcal{K}\left (\dfrac {r}{r'}\right )&\text{for }\quad r\lt r', \end{cases}\end{equation}

where $\mathcal{K}$ is the complete elliptic integral of the first kind. It is noted that as the argument approaches unity, the elliptic integral behaves like a logarithmic function, namely

(3.33) \begin{align} \mathcal{K}(r)\sim -\frac 12\ln (1-r)\quad \text{as } \quad r\rightarrow 1^-. \end{align}

Figure 3.

A true-colour image of the Indian Ocean taken by Terra’s moderate resolution imaging spectroradiometer captures arched cloud formations that indicate atmospheric gravity waves (Yiğit & Medvedev Reference Yiğit and Medvedev2019).

3.3. Canonical variables

It was proved by Benjamin & Bridges (Reference Benjamin and Bridges1997) that the two-layer interfacial wave problem is a Hamiltonian system, and interfacial variables $\xi ={\varPhi ^-}-{R\varPhi ^+}$ and $\eta$ form a canonical conjugate pair. Therefore, it is also common to write the governing equations in terms of canonical variables. First, it can be inferred from (3.9) that

(3.34) \begin{equation} \varPhi ^-=\xi +R\varPhi ^+=\xi -\mu \textit{RG}_0^{-1}\eta _t+O(\mu ^2). \end{equation}

Substituting (3.34) into (3.8), (3.11), and retaining terms valid up to $O(\mu )$ , one obtains

(3.35) \begin{equation} \left \{ \begin{aligned} &(\eta -b)_t+\boldsymbol{\nabla }\boldsymbol{\cdot }\big [(1+\eta -b)\,\boldsymbol{\nabla }\big (\xi -\mu \textit{RG}_0^{-1}\eta _t\big )\big ]=0,\\[3pt] &\xi _t+\eta +\frac 12|\boldsymbol{\nabla }\xi |^2-\mu R(\boldsymbol{\nabla }\xi )\boldsymbol{\cdot }\big (\boldsymbol{\nabla }G_0^{-1}\eta _t\big )=0, \end{aligned}\right . \end{equation}

which is a closed system for $\eta$ and $\xi$ . We can show that the above system is linearly stable, hence no further regularisation is needed. Actually, if we assume $b=0$ and $\eta ,\xi \ll 1$ , then the system can be linearised and simplified to

(3.36) \begin{equation} \eta _{tt}-\varDelta \eta +\mu \textit{RG}_0\eta _{tt}=0. \end{equation}

For a monochromatic wave $\eta =\widehat {\eta }\,\textrm{e}^{\textrm{i}(\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}-\omega t)}$ , the dispersion relation between wavenumber $\boldsymbol {k}$ and frequency $\omega$ is

(3.37) \begin{equation} \omega ^2=\frac {k^2}{1+\mu Rk}. \end{equation}

Recalling the full dispersion relation for two-layer interfacial waves with the semi-infinite upper layer, one can rewrite (2.15) in dimensionless form as

(3.38) \begin{eqnarray} \omega ^2=\frac {k\tanh (\mu k)}{\mu \left [1+R\tanh (\mu k)\right ]}. \end{eqnarray}

The right-hand side of (3.37) is an approximation of the full dispersion relation (3.38) for $\mu \ll 1$ . It is positive definite, and the phase speed $\omega /k$ is monotonically decreasing.

5.1. Linear theory

Before presenting numerical solutions for the nonlinear equations, we discuss the analytical results of linear theory. Linearising the system (4.16) and retaining terms valid up to $O(\mu )$ , one obtains

(5.1) \begin{equation} \eta _{tt}-\varDelta \eta +\mu \textit{RG}_0\eta _{tt}=b_{tt}, \end{equation}

where the two equations have been merged into a single one by eliminating $\xi$ . We assume that $b(x,y,t)$ is a Dirac delta function moving with constant speed $c$ in the $x$ -direction, namely $b=\delta (x+ct,y)$ , and look for the travelling-wave solution moving with the same speed. The problem can then be solved using the Fourier transform:

(5.2) \begin{equation} \eta =-\iint \frac {c^2k_x^2\,\textrm{e}^{\textrm{i}\boldsymbol {k}\boldsymbol{\cdot }\boldsymbol {r}}}{k^2-c^2k_x^2(1+\mu Rk)}\,\mathrm{d}\boldsymbol {k}. \end{equation}

We rewrite the above integral in the polar coordinate system using $\boldsymbol {k}=(k\cos \theta ,k\sin \theta )$ as

(5.3) \begin{equation} \eta =-\iint \frac {c^2k_x^2k\,\textrm{e}^{\textrm{i} xk(\cos \theta +\alpha \sin \theta )}}{k^2-c^2k^2\cos ^2\theta \,(1+\mu Rk)}\,\mathrm{d}k\,\mathrm{d}\theta, \end{equation}

where $\alpha =y/x$ represents the slope of a ray through the origin. The $k$ -integral is concentrated on the null points of $k^2-c^2k^2\cos ^2\theta (1+\mu Rk)$ , which gives

(5.4) \begin{equation} k=\frac {1}{\mu Rc^2}\left (\mathrm{sec}^2\theta -c^2\right ). \end{equation}

Therefore, for $\alpha$ fixed and $x$ large, the stationary phase gives the dominant contribution when $\varPhi (\theta ,c)=k(\theta ,c)\,(\cos \theta +\alpha \sin \theta )$ has a stationary point

(5.5) \begin{equation} \varPhi _\theta =0\;\Longrightarrow \;\alpha =-\frac {\mathrm{sec}\,\theta \tan \theta +c^2\sin \theta }{\mathrm{sec}\,\theta \,(2\tan ^2\theta +1)-c^2\cos \theta }, \end{equation}

where (5.4) has been used. For each ray $\alpha$ , there are waves with angle $\theta (k)$ . The extremum in $\alpha$ for which the relation (5.5) is true, denoted by $\alpha ^*$ , yields the ‘Kelvin’ wake angle $\arctan \alpha ^*$ . For example, $\alpha ^*=\arctan ({\sqrt {2}}/{4})$ for $c=0$ , and $\alpha ^*\rightarrow 90^\circ$ as $c\rightarrow 1$ . The curve of $\arctan \alpha ^*(c)$ is shown in figure 4(a) with a solid line. Accordingly, the critical wave angle and wavenumber for achieving the Kelvin wave angle are shown as a dashed curve in figure 4(a) and a solid line in figure 4(b), respectively, both of which are strictly decreasing functions of speed.

Figure 4.

(a) Kelvin wake angle $\arctan \alpha ^*$ versus speed (solid line), and critical $\theta$ for reaching the Kelvin wake angle (dashed line). (b) Critical wavenumber $k^*$ for reaching the Kelvin wake angle.

5.2. Wakes in the weakly nonlinear model

Weakly nonlinear long-wave theory requires a balance between weak dispersion and weak nonlinearity near the non-dispersive point ( $k=0$ ). The typical example is the KdV equation for pure gravity waves in shallow water. Along this line, we let $\epsilon =O (\mu )$ and ignore the terms of $O (\epsilon ^2,\mu ^2,\epsilon \mu )$ in (4.18). It is convenient to use $\xi$ as the unknown. Substituting (4.15) into (4.18) and retaining terms valid up to $O (\epsilon ,\mu )$ , one obtains

(5.6) \begin{equation} \begin{aligned} 0={}&\left (1+\mu \textit{RG}_0\right )\xi _{tt}-\varDelta \xi -\mu \textit{RG}_0^{-1}(\eta -b)_{ttt}-\mu \textit{RG}_0\eta _t+b_t\\ &{}+\epsilon \big [|\boldsymbol{\nabla }\xi |^2_t+\xi _t\,\varDelta \xi \big ]+\epsilon \boldsymbol{\nabla }\boldsymbol{\cdot }(b\,\boldsymbol{\nabla }\xi ). \end{aligned} \end{equation}

Upon noticing (4.16), we can infer that $-\mu \textit{RG}_0^{-1}(\eta -b)_{ttt}-\mu \textit{RG}_0\eta _t=O (\mu ^2,\epsilon \mu )$ , which is also negligible. Finally, returning to the original dependent and independent variables, (5.6) can be rewritten as

(5.7) \begin{align} \big [1+RhG_0\big ]\xi _{tt}-(1-R)gh\,\varDelta \xi +\big [|\boldsymbol{\nabla }\xi |^2_t+\xi _t\,\varDelta \xi \big]+(1-R)gb_t+(1-R)g\,\boldsymbol{\nabla }\boldsymbol{\cdot }(b\,\boldsymbol{\nabla }\xi )=0. \end{align}

We look for travelling-wave solutions to (5.6) or (5.7) by means of a numerical method to simulate the internal wave wake when wind passes over a local topographic high. The numerical scheme for constructing wake solutions is a generalisation of Petviashvili’s method (Petviashvili Reference Rottman and Einaudi1976), and the basic idea is to perform the fixed-point iteration in the Fourier space. Assume that the wind moves steadily with speed $c$ in the $x$ -direction, or equivalently, the bottom topography moves in the opposite direction with the same speed, namely $b(x,y,t)=b(x+ct,y)$ . Passing to the travelling coordinate system and applying the two-dimensional Fourier transform to (5.6), one obtains

(5.8) \begin{equation} \mathcal{F}[\xi]=-\frac {c\,\mathcal{F}[b_x]+\epsilon c\,\mathcal{F}\big [|\boldsymbol{\nabla }\xi |^2_x+\xi _x\,\varDelta \xi \big]+\epsilon\, \mathcal{F} [\boldsymbol{\nabla }\boldsymbol{\cdot }(b\,\boldsymbol{\nabla }\xi )]}{k^2-c^2k_x^2 (1+\mu Rk)}\triangleq \mathcal{P}[\mathcal{F}[\xi],\mathcal{F}[b]]. \end{equation}

Therefore, the fixed-point iteration process can be proposed as

(5.9) \begin{equation} \mathcal{F}\left [\xi _{n+1}\right ]=\mathcal{P}\left [\mathcal{F}\left [\xi _n\right ],\mathcal{F}\left [b\right ]\right ]. \end{equation}

In contrast to Petviashvili (Reference Rottman and Einaudi1976) and Ablowitz et al. (Reference Ablowitz, Fokas and Musslimani2006), no normalisation factor is needed due to the bottom topography, which is an external forcing to suppress the unlimited increase or decrease in amplitude during the iteration.

We then present numerical experiments of internal wave wakes. Wave profiles were computed in (5.8) with $R=0.99$ and the bottom topography

(5.10) \begin{equation} b(x,y)=\beta \,\textrm{e}^{-\gamma \left (x^2+y^2\right )}\,, \end{equation}

where $\beta$ and $\gamma$ are constants characterising its amplitude and width. A typical example is shown in figure 5(a) for a side view, figure 5(b) for an $x$ -cross-section, and figure 5(c) for a contour plot. The wave profile, which features a V-shaped pattern consisting of two wake lines, shows qualitative agreement with the field observation demonstrated in figure 1: each wake line is made up of ‘feathered’ waves of limited transverse length and at the same angle $\theta$ , and the interior of the wake is dominated by slightly curved wave fronts that are normal to the flow (see figure 5 c). Between the wake line and the interior is a transition region of shorter waves. Both linear and weakly nonlinear solutions are plotted together in figure 5(b) for comparison purposes, with the linear solution slightly overestimating the amplitude. It is noted that the linear solution was obtained by letting $\epsilon =0$ and taking the inverse Fourier transform of (5.8).

Figure 5.

The internal wave profile for $R=0.99$ , $\epsilon =\mu =0.2$ , $\beta =0.2$ , $\gamma =0.1$ , $c=0.91$ : (a) side view; (b) $x$ -cross-section at $y=0$ , and (c) contour plot. Only the active part is shown for better visibility in (b,c). The cross-section shown in (b) includes the solution to the nonlinear equation (blue solid curve), the solution to the linear equation (red dashed curve), and the bottom topography (black curve).

In computations, small artificial viscosity was introduced to break the symmetry of the solution about $x=0$ , and the wave profile is plotted by using the relation $\eta \approx -c\xi _x-({\epsilon }/{2})|\boldsymbol{\nabla }\xi |^2$ . Since wake solutions decay very slowly for large $c$ , a large number of mesh points are needed. For example, the numerical experiment shown in figure 5 was carried out for ${\rm d}x={\rm d}y=0.25$ and $2000\times 3200$ Fourier modes to obtain accurate enough solutions. This requirement reflects the advantage of the iterative scheme (5.8) over the traditional Newton–Raphson method. The Petviashvili iteration is fast and stable, avoids the construction of large Jacobian matrices, and incurs no additional computational cost when searching for non-symmetric steady-state solutions.

Figure 6.

Contour plots of internal waves wakes for incoming flow at various speeds and linear predictions of wake angles (dark dashed lines); the parameters are $R=0.99$ , $\epsilon =\mu =0.1$ , $\beta =0.1$ , $\gamma =4$ , with varying $c$ : (a) $c=0.85$ and $\alpha ^*=0.9649$ ; (b) $c=0.9$ and $\alpha ^*=1.2150$ ; (c) $c=0.95$ and $\alpha ^*=1.7697$ .

Contour plots of the internal wave wake for incoming flow at various velocities are shown in figures 6(a–c) for $c=0.85$ , $0.9$ and $0.95$ , respectively. The angle between wake lines increases with increasing wind speed toward the dimensionless critical value 1, akin to the ship wake on water of finite depth when gradually approaching the criticality (Havelock Reference Helfrich and Melville1908). As this speed is approached, ‘feathered’ waves on the wake line become more prominent, leading to the conjecture that the field observation in figure 1 has near-critical flow. The dashed lines are linear predictions (5.5), which show good agreement with the nonlinear solutions. Results of various aspect ratios of the bottom ridge are compared for a fixed upstream velocity ( $c=0.9$ , say). It is evident in figure 7(a) that increasing the ridge’s height leads to an amplification of wave amplitude. In addition, as the width enlarges (or equivalently, $\gamma$ decreases), the wave amplitude increases (see figure 7 b). As expected (and not shown), for a fixed speed $c$ , changing the vertical and horizontal sizes of the bottom topography has little impact on the wake angle.

Figure 7.

Effects of the aspect ratio of bottom topography on internal wave wakes obtained by simulating (5.8) with $R=0.99$ , $\epsilon =\mu =0.1$ and $c=0.9$ . (a) Green solid line indicates $\beta =0.05$ and $\gamma =4$ ; blue dashed line indicates $\beta =0.1$ and $\gamma =4$ ; red dotted line indicates $\beta =0.2$ and $\gamma =4$ . (b) Green solid line indicates $\beta =0.1$ and $\gamma =1$ ; blue dashed line indicates $\beta =0.1$ and $\gamma =4$ ; red dotted line indicates $\beta =0.1$ and $\gamma =16$ .

Figure 8.

An internal wave wake obtained by simulating (5.11) with $R=0.99$ , $\mu =0.2$ and $c=0.7$ : (a) side view and (b) top view.

5.3. Wakes in the strongly nonlinear model

It is also interesting to compute steady solutions in the strongly nonlinear equations developed in the present paper, which can be achieved by slightly generalising the Petviashvili method. Performing the two-dimensional Fourier transform on the system (3.8), (3.11) yields

(5.11) \begin{equation} \begin{pmatrix} \mathcal{F}[\eta ]\\[3pt] \mathcal{F}[\varPhi ^-] \end{pmatrix}= \begin{pmatrix} \text{i} ck_x & -k^2\\[3pt] 1-\dfrac {\mu Rc^2k_x^2}{k} & \text{i} ck_x \end{pmatrix}^{-1} \begin{pmatrix} \mathcal{F}\left [b_x-\boldsymbol{\nabla }\boldsymbol{\cdot }\left [(\eta -b)\,\boldsymbol{\nabla }\varPhi ^-\right ]\right ]\\[3pt] -\frac 12\mathcal{F}\left [|\boldsymbol{\nabla }\varPhi ^-|^2\right ] \end{pmatrix}\!. \end{equation}

An iterative algorithm can then be constructed as (5.9). Based on this algorithm, we first explored the relationship between the wave amplitude and the vertical size of the bottom topography. We chose $\gamma =16$ , $R=0.99$ , $\mu =0.2$ , $c=0.7$ , and considered $\beta$ as a bifurcation parameter. Nonlinear solutions were computed by varying the amplitude of the bottom topography. As a full exploration of the solution space lies beyond the present scope, the continuation process was terminated at $\beta \approx 0.39$ , where the fixed-point iteration failed to converge. This divergence likely signals a turning point in the solution branch – a regime where the Jacobian of the system becomes singular, rendering standard iterative schemes non-convergent. This limit primarily reflects a numerical constraint of the current solver; while arc length continuation could, in principle, bypass such singularities, its integration with the Petviashvili iteration is non-trivial. Further investigation using advanced continuation methods would be required to resolve the branch and explore potential secondary solutions beyond this threshold. The top and side views of the strongly nonlinear solution for $\beta =0.39$ are shown in figure 8, where the theoretical prediction of the wake angle (dashed line in figure 8 b) again agrees well with the numerical solution. For comparison purposes, in figure 9, we show the $x$ -cross-sections of wave profiles for the fully nonlinear model (blue solid curve) and the weakly nonlinear theory (red dashed curve) with different values of $\beta$ . In this regime, the weakly nonlinear solutions obtained from (5.8) always underestimate the wave amplitude, and as $\beta$ increases, a phase shift between the two solutions becomes apparent.

Figure 9.

Comparison of the $x$ -cross-section between the fully nonlinear solution (blue solid line) and the weakly nonlinear solution (red dashed line) for (a) $\beta =0.13$ , (b) $\beta =0.26$ , (c) $\beta =0.39$ .

Nonlinear effects are much stronger in flows past an isolated depression or crater, which can be easily computed up to large amplitudes. We show a strongly nonlinear numerical result far from criticality in figure 10. The theoretical prediction ( $\alpha ^*=0.4639$ for $c=0.5$ ) of the wake angle still makes good agreement (see figure 10 a). However, in figure 10(b), it is clear that the profiles are substantially different: the linear theory far overestimates the wave amplitude, and the far field is out of phase by approximately a half-wavelength.

Figure 10.

A nonlinear solution for subcritical flow past a submerged crater centred at the origin. The solution was simulated with $R=0.99$ , $\mu =0.2$ , $c=0.5$ , $\beta =-0.6$ and $\gamma =16$ . (a) Contour plot of the nonlinear solution and the linear prediction of the wake angle (dark dashed line with $\alpha ^*=0.4639$ ). (b) Comparison of the $x$ -cross-section between the fully nonlinear solution (blue solid line) and the linear prediction (red dashed line).

In certain regions, internal wave wakes generated by wind blowing over multiple islands may interact and superpose. A typical example is shown in figure 11(a), where a complex wake can be seen downwind of the South Sandwich Islands chain. Figure 11(b) shows a numerical result simulated with the bottom topography

(5.12) \begin{equation} \begin{aligned} b(x,y)={}&0.2\textrm{e}^{-16x^2-16y^2}+0.2\textrm{e}^{-16x^2-16(y-2.5)^2}+0.2\textrm{e}^{-16(x+1)^2-16(y+2.5)^2}\\ &{}+0.2\textrm{e}^{-16(x+1)^2-16(y-6)^2}+0.2\textrm{e}^{-16(x+3.5)^2-16(y+5)^2} \end{aligned} \end{equation}

to mimic the distribution of islands. A steady solution was obtained for $c=0.7$ , $R=0.99$ and $\mu =0.2$ , indicating that all islands trigger V-shaped cloud patterns that extend downwind, with wave–wave interactions occurring when two or more wake regions overlap. Qualitative aspects of the satellite image are reproduced despite our model being very simple, and in particular, assuming a (spatially) constant wind.

Figure 11.

(a) An interlocking series of mountain-wave clouds near the South Sandwich Islands in the southern Atlantic Ocean acquired by NASA’s Terra satellite on 5 February 2020. (b) A numerical simulation of the internal wave wake generated by flow over five separated uplifts given in (5.12), where the yellow dots represent centres of islands.

5.4. Ring waves

To simulate annular internal gravity waves as shown in figure 3, we focus on the weakly nonlinear equation (5.6) in the radially symmetric form

(5.13) \begin{equation} \left (1+\mu \textit{RG}_0\right )\xi _{tt}-\frac {1}{r}\left (r\xi _r\right )_r+\epsilon \left [\left (\xi _r^2\right )_t+\frac {\xi _t}{r}\left (r\xi _r\right )_r\right ]=0 \end{equation}

when ignoring the bottom topography. We use the pseudo-spectral method based on the Hankel transform to conduct the time-dependent numerical simulations of (5.13). The Hankel transform of order $\nu$ and its inverse are defined as

(5.14) \begin{align} &F(k)=\mathcal{H}_\nu [f(r)]=\int _0^\infty r\,f(r)\,J_\nu (kr)\mathrm{d}r,\nonumber \\ &f(r)=\mathcal{H}_\nu ^{-1}[F(k)]=\int _0^\infty k\,F(k)\,J_\nu (kr)\mathrm{d}k, \end{align}

respectively, where $J_\nu$ is the Bessel function of the first kind of order $\nu$ . The Hankel transform of order zero is equivalent to a two-dimensional Fourier transform of radial symmetry, which arises naturally in solving partial differential equations formulated in cylindrical coordinates. It follows from the definition of the Hankel transform of order zero and the definition of the Bessel function that

(5.15) \begin{equation} \begin{aligned} \left (\frac {\partial ^2}{\partial r^2}+\frac {1}{r}\frac {\partial }{\partial r}\right )f=&\int _0^\infty k\,\mathcal{H}_0[f](k)\left (\frac {\partial ^2}{\partial r^2}+\frac {1}{r}\frac {\partial }{\partial r}\right )J_0(kr)\,\mathrm{d}k\\ =&-\int _0^\infty k^3\,\mathcal{H}_0[f](k)\,J_0(kr)\mathrm{d}k. \end{aligned} \end{equation}

Equation (5.15) indicates that the symbol of the axisymmetric Laplace operator is $-k^2$ . We convert (5.13) into the spectral space, namely

(5.16) \begin{equation} \widehat {\xi }_{tt}+\omega ^2\widehat {\xi }=-\frac {\epsilon }{1+\mu Rk}\,\mathcal{H}_0\left [2\xi _r\left (\xi _t\right )_r+\frac {\xi _t}{r}\left (r\xi _r\right )_r\right ]\!, \end{equation}

where the hat denotes the Hankel transform of order zero, and $\omega ^2$ is given by (3.37). Following Milewski & Tabak (Reference Petviashvili1998), we can recast the system (5.13) into a first-order evolution equation by introducing $\widehat {p}=\widehat {\xi }_t-\text{i}\omega \widehat {\xi }$ for positive $\omega$ . Therefore, the system (5.16) can be rewritten as

(5.17) \begin{equation} \widehat {p}_t+\text{i}\omega \widehat {p}=-\frac {\epsilon }{1+\mu Rk}\,\mathcal{H}_0\left [2\xi _r\left (\xi _t\right )_r+\frac {\xi _t}{r}\left (r\xi _r\right )_r\right ]\!, \end{equation}

while $\widehat {\xi }_t$ and $\widehat {\xi }$ can be recovered by simply taking the real and imaginary parts of $\widehat {p}$ , respectively. The time integration of the complex equation (5.17) is implemented using integrating factors and the fourth-order Runge–Kutta method. The integral of the Hankel transform is numerically approximated by the combination of linear interpolation and direct integration. For a function with compact support or at least sufficiently fast decay, an equidistance radial grid $0=r_0\lt r_1\lt \cdots \lt r_{N-1}$ with sufficiently large $r_{N-1}$ can be introduced, and therefore

(5.18) \begin{equation} \begin{aligned} \mathcal{H}_0[f](k)\approx &\sum _{n=0}^{N-1}\frac {f(r_{n+1})+f(r_n)}{2}\int _{r_n}^{r_{n+1}}r\,J_0(kr)\,\mathrm{d}r\\[3pt] =&\sum _{n=0}^{N-1}\frac {f(r_{n+1})+f(r_n)}{2}\frac {r_{n+1}\,J_1(kr_{n+1})-r_n\,J_1(kr_n)}{k}. \end{aligned} \end{equation}

The nonlinear terms are computed in physical space, and spatial derivatives are approximated using the central difference scheme. Special treatments are required at the origin for both the Hankel transform and the Laplace operator, namely

(5.19) \begin{equation} \mathcal{H}_0[f](0)\approx \sum _{n=0}^{N-1}\frac {f(r_{n+1})+f(r_n)}{2}\frac {r_{n+1}^2-r_n^2}{2}\quad \text{and}\quad \lim _{r\rightarrow 0}\frac {1}{r}\frac {\partial f}{\partial r}=\frac {\partial ^2f}{\partial r^2}\bigg |_{r=0}. \end{equation}

The numerical algorithm was first validated by letting $\epsilon =\mu =0$ and comparing numerical results with the exact radial solution of the two-dimensional wave equation,

(5.20) \begin{equation} \xi (r,t)=\alpha \operatorname {Re}\left \{\frac {1+\text{i} t/\beta }{\left [(1+\text{i} t/\beta )^2+r^2/\beta ^2\right ]^{3/2}}\right \}\!, \end{equation}

where $\alpha$ and $\beta$ are arbitrary constants. This can be used to check the accuracy of the discrete Hankel transform. Further validation was carried out for non-zero $\epsilon$ and $\mu$ by comparing with the results obtained via simulating (5.6) with axisymmetric initial data (the interested reader is referred to Milewski & Tabak (Reference Petviashvili1998) for details of the numerical scheme). Good agreements were found in both types of comparisons.

Figure 12.

A radial dispersive shock wave generated by the initial conditions $\xi _t(r,0)=0$ and $\xi (r,0)=5/(1+r^2)^{3/2}$ . (a) Initial wave profile (blue solid line) and initial velocity potential (red dashed line). (b) Solution at $t=30$ : wave profile $\eta$ (blue solid line) and velocity potential $\xi$ (red dashed line). (c) The contour plot of the wave profile at $t=60$ .

The first numerical experiment was conducted on the concentric weakly nonlinear equation (5.13) with $\epsilon =\mu =0.1$ , $R=0.99$ , and the initial condition given by (5.20) for $\alpha =5$ and $\beta =1$ (see figure 12(a) for the initial $\eta$ and $\xi$ ). The wave profile at $t=30$ is shown in the blue solid curve in figure 12(b), featuring a well-developed internal dispersive wave train. In addition, the exact solution to the linear wave equation (black dot-dashed line) and the numerical solution to the weakly nonlinear weakly dispersive equation (red dashed line) for $\xi$ are compared in the same figure for $t=30$ . The numerical solution for $\xi$ also displays an oscillatory decaying wave train behind the leading trough. The top view of the radial dispersive wave train at $t=60$ is shown in the relief plot in figure 12(c), akin to the pattern formation of successive arched clouds seen in figure 3.

Figure 13.

Time evolution of a dam-break problem leads to two outward-propagating dispersive shock waves: (a) initial wave profile; (b) wave profile at $t=30$ ; (c) relief plot at $t=60$ .

Figure 14.

Two-way propagation of a radial internal wave: (a) initial wave profile; (b) wave profile at $t=20$ ; (c) relief plot at $t=50$ .

We finally consider the dam-break problem in the radially symmetric geometry. We consider a spatially confined pulse-type initial condition

(5.21) \begin{equation} \xi (r,0)=0,\quad \xi _t(r,0)=0.5\,\textrm{e}^{-36\left [(r-r_0)/l\right ]^{36}}, \end{equation}

where $r_0$ is the centre of the plateau, and $l$ is used to characterise its width. Numerical experiments were performed for $r_0=0$ and $r_0=100$ with $l=10$ , and the results are shown in figures 13 and 14, respectively. This can be considered as an approximation of a dispersive Riemann problem consisting of two annular constant states connected by an abrupt transition. For $r_0=0$ , two outward-propagating oscillatory ring structures corresponding to a dispersive shock and a modulated wave train are observed. A similar wave profile was also found in the concentric KdV equation for internal waves (Khusnutdinova & Zhang Reference Lindzen2016). However, the numerical result for $r_0=100$ shows different behaviour, with the appearance of two dispersive shock waves: one increasing in amplitude and propagating inwards (i.e. a dispersive shock tending to focus at the origin), and the other outward-propagating and decreasing in amplitude, as in the previous experiment. In addition, two modulated wave trains are generated between the two dispersive shocks. This structure arises from the bidirectional nature of (5.13), and cannot be captured by unidirectional models. Initially, the two middle wave trains are superposed, but eventually the structures separate as shown in figure 14(c).