2.1. Relationship between the fBL equation and the fKP equation

Consider a three-dimensional incompressible and inviscid fluid flowing irrotationally, bounded above by a free surface $z=\eta (x,y,t)$ and below by a varying bottom topography $b(x,y,t)$ , where $x$ and $y$ are the horizontal coordinates, and the $z$ -axis points upwards with $z=0$ the undisturbed free surface. The velocity potential on the free surface, $\xi (x,y,t)$ , is governed by the fBL equation (Milewski Reference Milewski1998)

(2.1) \begin{equation} \left (1-\frac {\mu ^2}{3}\unicode{x1D6E5} \right )\xi _{tt}-\unicode{x1D6E5} \xi +\epsilon \left [\left (|\nabla \xi |^2\right )_t+\xi _t\unicode{x1D6E5} \xi \right ]+\epsilon b_t=0\,. \end{equation}

Here, $\mu$ and $\epsilon$ are small parameters measuring the dispersive and nonlinear effects, respectively, and we assume $\epsilon =O(\mu ^2)$ to balance these two effects. We should point out that (2.1) is derived under the assumption of a small perturbation in bottom topography, specifically that the ratio of the bottom amplitude to the mean depth of the fluid is of order $\epsilon ^2$ . This equation is commonly used to explore the mechanisms by which waves are generated in flow over topography (Milewski & Tabak Reference Milewski and Tabak1999). It is also noted that (2.1) can also describe nonlinear internal waves, with variations primarily manifesting as differences in the coefficients of the constituent terms (Yuan et al. Reference Yuan, Wang and Chen2020).

Equation (2.1) can be reduced to the fKP equation by the following change of variables. Assuming a primarily left-propagating wave with weak transverse variation, we substitute

(2.2) \begin{equation} X=x+t\,,\quad Y=\mu y\,,\quad T=\epsilon t\,, \end{equation}

into (2.1) and neglect terms of order $\textit {O}(\epsilon ^2,\epsilon \mu )$ , leading to

(2.3) \begin{equation} \xi _{XT}+\frac {3}{2}\xi _{X}\xi _{XX}-\frac {\mu ^2}{\epsilon }\left (\frac {1}{6}\xi _{XXXX}+\frac {1}{2}\xi _{YY} \right )+\frac {1}{2}b_{X}=0\,. \end{equation}

Returning to the original independent variables and letting $\zeta =-\xi _x$ , (2.3) becomes the fKP equation,

(2.4) \begin{equation} \left (\zeta _t-\zeta _x-\frac {3\epsilon }{2}\zeta \zeta _x-\frac {\mu ^2}{6}\zeta _{xxx} \right )_x-\frac {1}{2}\zeta _{yy}=\frac {\epsilon }{2}b_{xx}\,, \end{equation}

which covers the KP equation under flat topography $b\equiv 0$ . Lee & Grimshaw (Reference Lee and Grimshaw1990) conducted a detailed investigation into the solution behaviour of (2.4) when a localized forcing disturbance moves with a near-critical speed in horizontally unbounded shallow water.

To normalize the fBL equation, we introduce new independent and dependent variables

(2.5) \begin{equation} \tilde {x}=\frac {\sqrt {3}}{\mu }x\,,\quad \tilde {y}=\frac {\sqrt {3}}{\mu }y\,,\quad \tilde {t}=\frac {\sqrt {3}}{\mu }t\,,\quad \tilde {\xi }=\frac {\sqrt {3}\epsilon }{\mu }\xi \,,\quad \tilde {b}=\frac {\epsilon \mu }{\sqrt {3}}b\,, \end{equation}

so that (2.1) can be transformed to the standard form

(2.6) \begin{equation} \xi _{tt}-\unicode{x1D6E5} \xi _{tt} -\unicode{x1D6E5} \xi +\left (|\nabla \xi |^2\right )_t+\xi _t\unicode{x1D6E5} \xi +b_t=0\,, \end{equation}

where all the tilde superscripts have been dropped for ease of notation. When the bottom topography is flat, the standard form of the fBL (2.6) reduces to the classic BL (1.6). In the subsequent sections, our analyses associated with the BL equation will be based on (2.6). Considering the flat topography situation, (2.6) admits the exact solution

(2.7) \begin{equation} \xi =\frac {a}{\delta \omega }\tanh \left [ \delta (kx+ly-\omega t) \right ]\,,\quad \delta =\frac {\sqrt {a}}{2|\omega |}\,, \end{equation}

where $k$ and $l$ denote the wavenumbers in the $x$ and $y$ directions, respectively, $\omega$ is the wave frequency and $a$ is a free parameter measuring the amplitude. These parameters are connected through the dispersion relation $a=\omega ^2/(k^2+l^2)-1$ . Correspondingly, under the scales (2.5), the fKP equation reads

(2.8) \begin{equation} \left (\zeta _t-\zeta _x-\frac {3}{2}\zeta \zeta _x-\frac {1}{2}\zeta _{xxx} \right )_x-\frac {1}{2}\zeta _{yy}=\frac {1}{2}b_{xx}\,. \end{equation}

For the flat bottom scenario, it admits the line soliton solution

(2.9) \begin{equation} \zeta ^{({KP})}=A^{({KP})}\mathrm{sech}^2\left [\frac {\sqrt {A^{({KP})}}}{2k} (kx+ly-\omega t) \right ]\,, \end{equation}

with the dispersion relation $A^{({KP})}=-2-l^2/k^2-2\omega /k$ . Finally, a simple rescaling

(2.10) \begin{equation} \tilde {t}=t\,,\quad \tilde {x}=-2^{\frac{1}{3}}x-2^{\frac{1}{3}}t\,,\quad \tilde {y}=2^{\frac{2}{3}}y\,,\quad \tilde {\zeta }=\frac {2^{\frac{2}{3}}}{3}\zeta \,, \end{equation}

shifts the fKP (2.8) with a flat bottom to the standard form (1.1), which serves as a bridge for the subsequent comparisons between the conclusions drawn from the KP equation, namely expressions (1.4) and (1.5), and the BL equation.

2.2. Numerical schemes

Here, we briefly describe the numerical methods for dealing with the fBL equation. The numerical algorithm is based on the pseudospectral method, and all the derivatives and pseudodifferential operators are computed with the aid of their Fourier symbols. The time integration of the equations is carried out using the classic fourth-order Runge–Kutta scheme with the method of integrating factors (see Milewski & Tabak (Reference Milewski and Tabak1999) for details). Firstly, to fit these waves in a periodic domain of computation, we need to make $\xi$ periodic, which can be achieved by writing

(2.11) \begin{equation} \Psi =\xi (\boldsymbol{x},t)+\boldsymbol{C}\cdot \boldsymbol{x}\,, \end{equation}

where the constant vector $\boldsymbol{C}=(C_1,C_2)$ is chosen to ensure that $\Psi$ is periodic. Thus, the fBL (2.6) can be rewritten as

(2.12) \begin{equation} \left ( 1-\unicode{x1D6E5} \right )\Psi _{tt} -2\boldsymbol{C}\cdot \nabla \Psi _t-\unicode{x1D6E5} \Psi =\mathcal{N}(\Psi, b)\,, \end{equation}

where $\mathcal{N}(\Psi, b)$ includes all the nonlinear and topography terms, shown as

(2.13) \begin{equation} \mathcal{N}(\Psi, b)=-|\nabla \Psi |^2_t-\Psi _t\unicode{x1D6E5} \Psi -b_t\,. \end{equation}

Secondly, by performing a Fourier transform of (2.12) with respect to the spatial variables, we can reformulate it into a first-order single complex evolution equation through the introduction of

(2.14) \begin{equation} \widehat {v_1}=\left ( \partial _t+\mathcal{S}+\mathrm{i} \mathcal{L} \right )\widehat {\Psi }\,,\quad \widehat {v_2}=\left ( \partial _t+\mathcal{S}-\mathrm{i} \mathcal{L} \right )\widehat {\Psi }\,, \end{equation}

with

(2.15) \begin{equation} \mathcal{S}=-\frac {\boldsymbol{C}\cdot \mathrm{i}\boldsymbol{k}}{1+|\boldsymbol{k}|^{2}}\,,\quad \mathcal{L}=\sqrt {\frac {|\boldsymbol{k}|^{2}}{1+|\boldsymbol{k}|^{2}}+\frac {\left (\boldsymbol{C}\cdot \boldsymbol{k}\right )^2}{\left (1+|\boldsymbol{k}|^{2}\right )^{2}}}\,, \end{equation}

where the hat denotes Fourier-transformed variables and $\boldsymbol{k}= ( k, l )$ is the wavenumber vector. Therefore, (2.12) can be recast as

(2.16) \begin{equation} \left (\partial _t+\mathcal{S}-\mathrm{i}\mathcal{L}\right )\widehat {v_1}=\frac {1}{1+|\boldsymbol{k}|^{2}}\widehat {\mathcal{N}}\,,\quad \left (\partial _t+\mathcal{S}+\mathrm{i}\mathcal{L}\right )\widehat {v_2}=\frac {1}{1+|\boldsymbol{k}|^{2}}\widehat {\mathcal{N}}\,. \end{equation}

Upon noticing that $\Psi$ is real so that $\widehat {\Psi }^*(\boldsymbol{k})=\widehat {\Psi }(-\boldsymbol{k})$ , the two equations in (2.16) are equivalent; thus, $\widehat {\Psi }$ and $\widehat {\Psi _t}$ can be recovered from $\widehat {v_1}$ alone with

(2.17a) \begin{align} &\qquad \widehat {\Psi }(\boldsymbol{k})=\frac {1}{2\mathrm{i}\mathcal{L}}\left [\widehat {v_1}(\boldsymbol{k})-\widehat {v_1}^*(-\boldsymbol{k})\right ]\,, \end{align}

(2.17b) \begin{align} &\widehat {\Psi }_t(\boldsymbol{k})=\frac {1}{2}\left [ \widehat {v_1}(\boldsymbol{k})+\widehat {v_1}^*(-\boldsymbol{k}) \right ]-\mathcal{S}(\boldsymbol{k})\widehat {\Psi }(\boldsymbol{k})\,, \end{align}

where the asterisk represents the complex conjugation. When $\boldsymbol{k}=(0,0)$ , it is noted that $\mathcal{L}(\boldsymbol{0})=0$ , which leads to the singularity of $\widehat {\Psi }(\boldsymbol{k})$ in (2.17a ). Indeed, the value of $\widehat {\Psi }(\boldsymbol{0})$ can be obtained by directly solving (2.17b ). The approach of integrating factors, in conjunction with the fourth-order Runge–Kutta scheme and the pseudospectral method, is utilised for the time stepping of the first equation in (2.16). Interested readers are referred to Milewski & Tabak (Reference Milewski and Tabak1999) for more details.

To investigate the non-periodic and non-decaying solutions to the fBL equation, a windowing method, as proposed by Schlatter et al. (Reference Schlatter, Adams and Kleiser2005), is implemented to facilitate successful numerical computations. Indeed, this method has been widely employed to study nonlinear interactions among solitary waves by Kao & Kodama (Reference Kao and Kodama2012), Ablowitz & Curtis (Reference Ablowitz and Curtis2013), Yuan et al. (Reference Yuan, Grimshaw, Johnson and Chen2018) and Yuan & Wang (Reference Yuan and Wang2022). For line solitons, we extend the computational domain sufficiently long in the $x$ direction to ensure that all solitons remain spatially distant from the boundary lines. Consequently, this expansion necessitates that these waves intersect with the boundary lines in the $y$ direction. To retain the effectiveness of the Fourier transform in the $y$ direction, a fundamental approach involves introducing a smooth function within the windowing method. This function exhibits compact support, with values that divide the solution into the interior and boundary layer regions. The smooth function $W(y)$ can be simply defined as

(2.18) \begin{equation} W(y)=\exp \left (-a\Big |\frac {y}{L_y}\Big |^n\right )\,, \end{equation}

where $y=\pm L_y$ denotes the boundary lines in the $y$ direction with the parameters $n=27$ and $a=(1.111)^n\ln 10$ throughout this paper (see Kao & Kodama (Reference Kao and Kodama2012) for more details). Correspondingly, at each time step, the solution within the computational domain is redefined as

(2.19) \begin{equation} \Psi =W\Psi +(1-W)\widetilde {\Psi }\,, \end{equation}

where $\Psi$ is the solution to the fBL equation and $\widetilde {\Psi }$ is assumed to satisfy the fBL equation in the boundary layers. Applying the windowing method obviates the requirement for the constant vector $\boldsymbol{C}$ to guarantee the periodicity of the boundaries in the $y$ direction. Thus, we take $C_2=0$ in (2.12). Submitting the transformation (2.19) into the fBL (2.12) and denoting $\psi (x,y,t)=W(y)\Psi (x,y,t)$ , an alternative equation for $\psi$ can be derived as follows:

(2.20) \begin{equation} (1-\unicode{x1D6E5} )\psi _{tt}-2C_1\psi _{xt}-\unicode{x1D6E5} \psi =\mathcal{N}(\psi, b)+G_1(W,\widetilde {\Psi },b)+G_{2}(W,\psi, \widetilde {\Psi },b)\,, \end{equation}

where

(2.21) \begin{align} G_1(W,\widetilde {\Psi },b)=&-\left (1-W\right )\left [ (1-\unicode{x1D6E5} )\widetilde {\Psi }_{tt}-2C_1\widetilde {\Psi }_{xt}-\unicode{x1D6E5} \widetilde {\Psi }-\mathcal{N}(\widetilde {\Psi },b)\right ]\,,\nonumber \\ G_{2}(W,\psi, \widetilde {\Psi },b)=&-\left ( W_{yy}\widetilde {\Psi }_{tt}+2W_y\widetilde {\Psi }_{ytt}+W_{yy}\widetilde {\Psi }+2W_y\widetilde {\Psi }_y \right )-(1-W)\left [ 2\left ( \nabla \psi \cdot \nabla \widetilde {\Psi } \right )_t \hphantom{\left ( |\nabla \widetilde {\Psi }|^2_t+\widetilde {\Psi }_t\unicode{x1D6E5} \widetilde {\Psi } \right )\widetilde {\Psi }-W_{yy}\widetilde {\Psi }_t \widetilde {\Psi }}\right .\nonumber \\ &-W\!\left .\left ( |\nabla \widetilde {\Psi }|^2_t+\widetilde {\Psi }_t\unicode{x1D6E5} \widetilde {\Psi } \right ) -2W_y\!\left ( \widetilde {\Psi }_y \widetilde {\Psi }_t+\widetilde {\Psi } \widetilde {\Psi }_{yt} \right )+\psi _t\unicode{x1D6E5} \widetilde {\Psi }\!-W_{yy}\widetilde {\Psi }_t \widetilde {\Psi } \right ]\nonumber \\ &-2W_y^2\widetilde {\Psi } \widetilde {\Psi }_t+ \! 2W_y \!\left ( \psi _y \widetilde {\Psi }_t+\psi _{yt}\widetilde {\Psi }+\psi _t \widetilde {\Psi }_y \right )+W_{yy}\psi _t \widetilde {\Psi }+(1-W)b_t\,, \end{align}

subject to the initial condition

(2.22) \begin{align} \psi (x,y,0)=\Psi (x,y,0)-(1-W)\widetilde {\Psi }(x,y,0)\,. \end{align}

Given the assumption that $\widehat {\Psi }$ complies with the fBL (2.12) within the boundary layers, we can deduce that $G_1=0$ , leaving only the term $G_2$ to be influential. To ensure the spectral accuracy of the results, we analytically calculate the terms related to $\widehat {\Psi }$ and $W$ in $G_2$ .

Based on the above numerical method, computations for (2.6) are performed with a computational domain of $x\times y=[-500,500]\times [-1000,1000]$ , discretized with $1024\times 2048$ Fourier modes along the propagating and transverse directions, respectively. The pseudospectral scheme is coupled with the 2/3 rule dealiasing technique (i.e. using a buffer of half the number of modes in each horizontal direction for each time step). The numerical results for the BL equation will be presented in a reference frame that moves in the $x$ direction at a speed $v$ . In the context of solitary wave interactions, the frame speed is chosen so that a single oblique solitary wave remains stationary when there is no wave interaction.

4.1. Mach expansion

The occurrence of Mach expansion originates from the bent soliton initial condition (4.1) with $\theta _0\lt 0$ , depicted in figure 2(a), and the negative sign of $\theta$ yields two rarefaction waves emanating from the point $(0,0)$ . Figure 2(a–c) show the evolution of the initial data (4.1) with $\omega _0=-1.1$ and $\theta _0=-15^\circ$ . In figures 2(b) and 2(c), a vertical soliton is observed in conjunction with the evolution of the bent soliton, identified as the Mach stem.

The specific characteristics of the Mach stem can be described by the hyperbolic system (4.2) associated with the initial data (4.1). In this case, one wave propagates along the positive direction of the $y$ axis with the Riemann variable $r$ held constant, while the other wave propagates along the negative direction of the $y$ axis with the other Riemann variable $s$ held constant. Therefore, the slope $q_w$ and the amplitude $a_w$ of the Mach stem are subject to the constraint

(4.7) \begin{equation} r(a_0,\tan \theta _0)=r(a_w,q_w)\,,\quad s(a_0,-\tan \theta _0)=s(a_w,q_w)\,, \end{equation}

which yields the expressions

(4.8) \begin{equation} q_w=0\,,\quad \theta _0=I_{a_0}^{a_w}=\int _{a_0}^{a_w} \frac {\sqrt {12a^2+25a+15}}{2(1+a)\sqrt {a}\sqrt {5+3a}}\text{d}a\,. \end{equation}

Here, the amplitude $a_w$ of the Mach stem can be obtained by numerically solving the following boundary value problem:

(4.9) \begin{equation} \frac {\text{d}a(\theta )}{\text{d}\theta }=-\frac {2(1+a)\sqrt {a}\sqrt {5+3a}}{\sqrt {12a^2+25a+15}}\,,\quad a(\theta _0)=a_0\,,\quad a(0)=a_w\,. \end{equation}

Figures 3(a) and 3(b) show the evolution of the amplitude $a$ and the slope $q$ for the initial data (4.1) with $\theta _0=-15^\circ$ and $a_0=0.21$ . Results depicted by the dark solid curves are derived from numerical simulations for the BL equation (1.6), whereas the outcomes illustrated by the red dashed curves are obtained from numerical simulations for the modulation equations (4.2). It can be observed from figures 3(a) and 3(b) that the amplitude and the slope are well predicted by the solutions of the modulation equations (4.2). Furthermore, after a long period of evolution, the amplitude and the slope of the Mach stem converge to $a_w\approx 0.088$ and $q_w=0$ , respectively, as derived from (4.8). Interestingly, we have observed that the Mach wave takes the form of a soliton, accompanied by a series of arc-shaped oscillations trailing behind the wavefront, which will be explained in § 6. The waveform of the Mach stem can be expressed as

(4.10) \begin{equation} a=a_w \text{sech}^2\left [ \frac {\sqrt {a_w}}{2|\omega _0|}(x-\omega _0t) \right ]\,,\quad \omega _0=-\sqrt {a_0+1}\,. \end{equation}

Furthermore, the simple wave solution to the modulation (4.2) with the bent soliton initial data (4.1) reveals

(4.11) \begin{equation} \begin{aligned} I_0^a&=\left \{ \begin{array}{ll} I_0^{a_0}\,, & |y|\gt V_0t\,,\\[2pt] \frac {\theta _0}{V_w-V_0}\frac {|y|}{t}-\frac {V_0I_0^{a_w}-V_wI_0^{a_0}}{V_w-V_0}\,, & V_wt\le |y|\le V_0 t\,,\\[8pt] I_0^{a_w}\,, & |y|\lt V_wt\,, \end{array}\right .\\[3pt] q&=\text{sgn}(y)\left \{ \begin{array}{ll} \tan \theta _0\,, & |y| \gt V_0 t\,,\\[3pt] \tan (\theta _0+I_{a}^{a_0})\,, & V_w t\le |y|\le V_0 t\,,\\[3pt] 0\,, & |y|\lt V_wt\,, \end{array}\right . \end{aligned} \end{equation}

with

(4.12) \begin{equation} V_w=V(a_w,0)\,,\quad V_0=V(a_0,\tan \theta _0)\,. \end{equation}

The amplitude of the Mach stem, as determined by the modulation (4.2), is illustrated in figure 3(c) through the contour plot. In the plot, the black dash–dotted lines and blue dashed lines represent the edges of the rarefaction wave, with the slopes given by $\pm V_0$ and $\pm V_w$ .

The appearance of the Mach stem is due to $a_w\gt 0$ . When $a_w=0$ , we can determine the critical angle of the slope for the occurrence of the Mach stem, which is given by

(4.13) \begin{equation} \theta _{{cr}}=\int _{a_0}^{0}\frac {\sqrt {12a^2+25a+15}}{2(1+a)\sqrt {a}\sqrt {5+3a}}\;\text{d}a\,. \end{equation}

Thus, in the case of $\theta _{{cr}}\le \theta _0\le 0$ , the bent soliton will evolve into a Mach stem perpendicular to the propagation direction, exhibiting a phenomenon similar to Mach expansion in gas dynamics.

Figure 3.

(a–b) Evolutions of the amplitude $a$ and slope $q$ for the bent soliton initial data (4.1) with $\theta _0=-15^\circ$ and $a_0=0.21$ at different times. They are obtained from the numerical simulations for the BL equation (1.6) (dark solid curves) and the modulation equations (4.2) (red dashed curves). (c) Space–time contour plots of the amplitude of the Mach stem with the same initial data, where the blue dashed lines and black dash–dotted lines characterise the edges of the evolution, predicted by the solutions (4.11).

4.2. Regular expansion

When $\theta _0\lt \theta _{{cr}}$ (corresponding to the negative value of $a_w$ ), the Mach stem disappears; instead, the regular expansion phenomenon occurs. Similar to the analysis in Ryskamp et al. (Reference Ryskamp, Maiden, Biondini and Hoefer2021), we let $a_w=0$ , resulting in $q_w\neq 0$ . By symmetry, we can obtain $q$ locally based on the simple wave criterion,

(4.14) \begin{equation} q_{w+}=\lim _{y\rightarrow 0^+}q(y,t)=-\lim _{y\rightarrow 0^-}q(y,t)=q_{w-}\,, \end{equation}

with the initial discontinuity at $y=0$ . Then, with the aid of the Riemann invariants, the solution for the rarefaction wave is given by

(4.15) \begin{equation} \begin{aligned} I_0^a&=\left \{ \begin{array}{ll} I_0^{a_0}\,, & |y| \gt V_0t\,,\\[3pt] \frac {I_0^{a_0}}{V_0-V_w}\frac {|y|}{t}-\frac {V_wI_0^{a_0}}{V_0-V_w}\,, & V_w t\le |y|\le V_0t\,,\\[6pt] 0\,, & |y|\lt V_w t\,, \end{array} \right .\\[2pt] q&=\text{sgn}(y)\left \{ \begin{array}{ll} \tan \theta _0\,, & |y|\gt V_0t\,,\\[2pt] \tan (\theta _0+I_{a}^{a_0}) & V_wt\le |y|\le V_0t\,,\\[2pt] q_{w+}\,, & |y|\lt V_wt\,, \end{array}\right . \end{aligned} \end{equation}

with

(4.16) \begin{equation} V_w=V(0,q_{w+})\,,\quad V_0=V(a_0,\tan \theta _0)\,. \end{equation}

Figure 4.

Time evolution of the bent-stem soliton initial data (4.1) with $\theta _0=-60^\circ$ and $a_0=0.21$ : (a) $t=0$ , (b) $t=500$ and (c) $t=1000$ . The snapshots are shown in a frame of reference moving with $v=-1.1/ \cos (-60^\circ )=-2.2$ in the $x$ direction.

Figure 4 shows the regular reflection of the bent soliton with $\theta =-60^\circ$ and $a_0=0.21$ for the initial condition (4.1). To mitigate the interference of trailing waves from the cyclic boundary on the outcomes, we extend the computational domain to $x\times y=[-2000,2000]\times [-1000,1000]$ in this scenario while maintaining Fourier modes unchanged. Unlike the Mach expansion shown in figure 2, a stable Mach stem formation does not emerge in this situation. What is noted instead is the presence of an arc-shaped wavefront, which exhibits a gradual reduction in amplitude attenuating with time and serves to connect two bent solitons.

4.3. Comparison between the BL and KP equations

To capture the effects of isotropism, we investigate wave propagation in both the BL and KP equations, and the main results are shown in figure 5, where we take the bent soliton in the form of (4.1) as the initial data. In fact, due to its reduction of the BL equation, the KP equation shares characteristics similar to those of the BL equation to some extent, especially for the central part of the soliton evolution (Ablowitz & Curtis Reference Ablowitz and Curtis2011; Gidel et al. Reference Gidel, Bokhove and Kalogirou2017; Yuan & Wang Reference Yuan and Wang2022). Thus, we mainly focus on the differences between these two equations for the Mach expansion of the oblique soliton, namely, the wall reduction factor and the speed of the Mach stem.

According to the modulation equations (3.1), the corresponding amplitude and critical slope of the Mach stem in the KP equation can be obtained as follows:

(4.17) \begin{equation} A^{({KP})}_w=\left ( \sqrt {A^{({KP})}_0}+\frac {q_0}{\sqrt {3}}\right )^2\,,\quad q^{({KP})}_{{cr}}=-\sqrt {3A^{({KP})}_0}\,, \end{equation}

where the initial amplitude $A^{({KP})}_0$ is derived based on the dispersion relation

(4.18) \begin{equation} A^{({KP})}_0=-2-q_0^2-2\omega _0/k=-2-\tan ^2 \theta _0-2\frac {\omega _0}{\cos \theta _0}\,. \end{equation}

Thus, one obtains the wall reduction factor, speed and critical angle with respect to the initial angel $\theta _0$ and speed $\omega _0$ , shown as

(4.19) \begin{equation} \begin{aligned} &\lambda ^{({KP})}=\frac {A^{({KP})}_w}{A^{({KP})}_0}=\left (\frac {\sqrt {2}\sin \theta _0}{\sqrt {3}\sqrt {-3-\cos 2\theta _0-4\omega _0\cos \theta _0}}+1\right )^2\,,\\ &c_w^{({KP})}=-\frac {A^{({KP})}_w}{2}-1\,,\\ &\theta _{{cr}}^{({KP})}=\arctan q^{({KP})}_{{cr}}=-\text{arcsec}\left (\frac {\sqrt {9\omega _0^2-8}-3\omega _0}{4}\right )\,, \end{aligned} \end{equation}

where the wall reduction factor refers to the ratio of the amplitude of the Mach stem to that of the incident soliton (see, Ryskamp et al. Reference Ryskamp, Maiden, Biondini and Hoefer2021).

To facilitate a comparative analysis between the BL and KP equations, based on the relationship $\zeta =-\xi _x$ , we define the soliton amplitude in the BL equation as $A^{({BL})}=-ka/\omega$ , with the dispersion relation $a=\omega ^2/( k^2+l^2 )-1$ . Thus, we can reveal the wall reduction factor, speed and critical angle for the BL equation,

(4.20) \begin{align} \lambda ^{({BL})}& =\frac {A^{({BL})}_w}{A^{({BL})}_0}=\frac {a_w\sqrt {a_0+1}}{k_0a_0\sqrt {a_w+1}}\,, \nonumber \\ c_w^{({BL})} & =-\sqrt {A^{({BL})}_w+1}\,,\nonumber\\ \theta _{{cr}}^{({BL})}& =\int _{a_0}^{0} \frac {\sqrt {12a^2+25a+15}}{2(1+a)\sqrt {a}\sqrt {5+3a}}\,\text{d}a\sim \sqrt {6\varepsilon }\left ( 1-\frac {11}{180}\varepsilon +\frac {2123}{36000}\varepsilon ^2 \right )+O(\varepsilon ^3)\,, \end{align}

with $a_0=\omega _0^2-1$ , $\varepsilon =-\omega _0-1$ and $a_w$ being determined through numerical solutions of (4.9).

Figure 5.

Comparative plot of the wall reduction factor (a) and the speed (b) of the Mach stem in the KP (red dashed curves) and BL (black solid curves) equations with the initial angle $\theta _0$ and speed $\omega _0$ of soliton. The dash–dotted curves denote the critical angle of the Mach expansion versus the initial speed $\omega _0$ for the KP (red) and BL (black) equations given by (4.19) and (4.20), respectively.

Figure 5 presents a comparison of the wall reduction factor and the speed of the Mach stem between the KP and BL equations as functions of the initial angle $\theta _0$ and speed $\omega _0$ . The wall reduction factor and the Mach stem speed for the KP and BL equations are depicted by the red dashed and black solid curves, respectively. The dash–dotted curves correspond to the critical angle of the Mach wave fan for both the KP and BL equations, derived from expressions (4.19) and (4.20). Clearly, compared with the KP equation, the decrease of the critical angle for the BL equation is more rapid as $|\omega _0|$ is far away from the unity. Upon visual inspection of figure 5, it is discernible that when the speed $|\omega _0|$ is close to unity, there is a pronounced correspondence between the two datasets. Furthermore, as $|\omega _0|$ increases, the divergence between the datasets becomes progressively more pronounced. Examination of figure 5(a) reveals that, for each speed, the wall reduction factor diminishes as the initial angle approaches the critical angle. This signifies that the amplitude of the Mach stem formed by the Mach expansion relative to the incident wave amplitude decreases. Additionally, it is observed that the wall reduction factor for the KP equation is consistently lower than that for the BL equation. This discrepancy becomes increasingly evident as the initial speed increases. While in figure 5(b), for the speed of the Mach stem $c_w$ , regardless of the initial speed, the differences in $c_w$ between the two models are insignificant as the initial angle surpasses the criticality. However, as $\theta _0$ approaches the critical angle $\theta _{{cr}}$ , a distinct divergence emerges in the values of $c_w$ .

5.1. Mach reflection

When $\theta _0\gt 0$ in the initial data for the bent soliton, not only does the Mach stem occur, but two solitons also arise symmetrically with the slope sgn $(y)q\lt 0$ . Figure 6 illustrates the evolution of the reverse bent soliton with $a_0=0.21$ and $\theta _0=20^\circ$ in the initial condition (4.1). In figure 6, one can observe that, as the reverse bent soliton $(a_0,q_0)$ propagates, the Mach stem $(a_w,q_w)$ and the soliton $(a_i,q_i)$ are generated at the bend point, and as time evolves, these three solitons maintain their intersection at the same point $P$ , which is identified as the Miles resonant solitons (Miles Reference Miles1977b ).

In this case, two oblique solitons (analogous to oblique shock waves in gas dynamics) symmetrically emanate from the bend point. Ryskamp et al. (Reference Ryskamp, Hoefer and Biondini2022) investigated a similar initial condition for the KP equation. Building upon the methodology outlined in Ryskamp et al. (Reference Ryskamp, Hoefer and Biondini2022), we employ the modified RH conditions to scrutinize the dynamics of the reverse bent soliton, shown as

(5.1) \begin{equation} s_y [\![\!\kern2pt f \kern2pt \!]\!]_{m,n} = [\![h\!]\!] _{m,n}\,, \end{equation}

where $s_y$ denotes the velocity of the discontinuity in the $y$ direction, the multivalued functions $f$ and $h$ satisfy the conservation law $f_t+h_y=0$ and, $m$ as well as $n$ , denote the count of values of the multivalued functions on either side of the discontinuity point.

Figure 6.

Snapshot during the evolution of the bent-stem soliton initial data (4.1) with $\theta _0=20^\circ$ and $a_0=0.21$ : (a) $t=0$ , (b) $t=300$ and (c) $t=600$ . The snapshots are shown in a frame of reference moving with $v=-1.1/\cos (20^\circ )\approx -1.1706$ in the $x$ direction.

Therefore, the modulation equation (3.10d ) can be represented in the conservation form as demonstrated by

(5.2) \begin{equation} f_t+g_x+h_y=0\,, \end{equation}

with

(5.3) \begin{equation} f=a\sqrt {a}\frac {5+4a}{1+a}\,,\quad g=-a\sqrt {a}\frac {5+3a}{\sqrt {1+a}}\frac {1}{\sqrt {1+q^2}}\,,\quad h=-a\sqrt {a}\frac {5+3a}{\sqrt {1+a}}\frac {q}{\sqrt {1+q^2}}\,. \end{equation}

Thus, the velocity $(s_x,s_y)$ at the discontinuity point $P$ satisfies the modified RH conditions, shown as

(5.4) \begin{equation} s_x=\frac {g_0+g_i-g_w}{f_0+f_i-f_w}\,,\quad s_y=\frac {h_0+h_i-h_w}{f_0+f_i-f_w}\,, \end{equation}

where the subscripts $w$ , $i$ and 0 represent the parameters of the Mach, reflected and incident solitons, respectively, as shown in figure 6(c)

It is noted that the KP equation is an integrable system with an elegant mathematical structure and infinite conservation laws. However, the BL equation is non-integrable; thus, some conclusions in Ryskamp et al. (Reference Ryskamp, Hoefer and Biondini2022) drawn from the KP equation may not directly apply to the BL equation. The following modifications are applied to address the difficulty of non-integrability. Considering the consistency equation (3.10c ), by alternately neglecting the $x$ and $y$ directions, we can decompose it into two separate equations in the conservation form

(5.5) \begin{equation} F^{(1)}_t+G^{(1)}_x=0\,,\quad F^{(2)}_t+H^{(2)}_y=0\,, \end{equation}

with

(5.6) \begin{equation} \begin{aligned} F^{(1)}&=\frac {1}{q}\,,\quad G^{(1)}=-\frac {\sqrt {1+q^2}\sqrt {1+a}}{q}\,,\\ F^{(2)}&=q\,,\quad H^{(2)}=-\sqrt {1+q^2}\sqrt {1+a}\,. \end{aligned} \end{equation}

Moreover, the velocity $(s_x,s_y)$ at the point $P$ satisfies the classical RH conditions

(5.7) \begin{equation} \begin{aligned} s_x=\frac {G^{(1)}_0-G^{(1)}_w}{F^{(1)}_0-F^{(1)}_w}=\frac {G^{(1)}_i-G^{(1)}_w}{F^{(1)}_i-F^{(1)}_w}=\frac {G^{(1)}_0-G^{(1)}_i}{F^{(1)}_0-F^{(1)}_i}\,,\\ s_y=\frac {H^{(2)}_0-H^{(2)}_w}{F^{(2)}_0-F^{(2)}_w}=\frac {H^{(2)}_i-H^{(2)}_w}{F^{(2)}_i-F^{(2)}_w}=\frac {H^{(2)}_0-H^{(2)}_i}{F^{(2)}_0-F^{(2)}_i}\,, \end{aligned} \end{equation}

where the fourth term of the continued equality in each expression can be derived from the second and third terms.

Figure 7.

(a) Evolutions of the amplitude $a$ for the bent soliton initial data (4.1) with $\theta _0=20^\circ$ and $a_0=0.21$ at $t=600$ . They are obtained from the direct numerical simulations of the BL equation (1.6) (dark solid curves) and the theoretical prediction based on the algebraic equations (5.4) and (5.7) (red dashed curves) with $a_w\approx 0.478149$ and $a_i\approx 0.0568022$ . (b) The corresponding contour plot of (a) (see also figure 6 c): the black dash–dotted lines represent the theoretical prediction obtained through the algebraic equations (5.4) and (5.7) with $q_i\approx 0.797368$ .

By solving the system of equations (5.4) and (5.7) and considering the initial data $(a_0,q_0)$ , we can derive the numerical solution $(a_w,q_w)$ and $(a_i,q_i)$ . Additionally, the velocity $(s_x,s_y)$ can also be accurately determined. The comparative illustration between the results of theoretical predictions and numerical simulations is depicted in figure 7. In figure 7(a), the dark solid lines correspond to numerical simulations of the BL equation (1.6), while the red dashed lines represent theoretical predictions derived from the algebraic expressions (5.4) and (5.7). The numerically determined values for the amplitude $a$ and the slope $q$ exhibit a remarkable concordance with the theoretical expectations. The contour plot, as depicted in figure 7(b), further substantiates the precision of the theoretical forecasts; the black dash–dotted lines therein correspond to the theoretical predictions culled from the aforementioned algebraic relations. The concurrence between the outcomes of the numerical simulations and the theoretical prognoses is striking, with the modulation theory adeptly delineating the evolution of reverse bent solitons, known as Mach reflection.

5.2. Regular reflection

When the velocity $s_y\leqslant 0$ , the Mach stem will not appear. In this case, by substituting $s_y=0$ into (5.4) and (5.7), the critical slope $q_{{cr}}$ is attained, which only depends on the initial amplitude $a_0$ . The corresponding parameters for the reflected soliton are as follows:

(5.8) \begin{equation} a_i=a_0\,,\quad q_i=-q_0\,. \end{equation}

Moreover, based on the expressions of $s_x$ in (5.4) and (5.7), one obtains the relation between $a_w$ and $a_0$ as

(5.9) \begin{equation} 2\left (4a_0+5\right )a_0^{\frac{3}{2}}a_w-\left (a_0+1\right )a_w^{\frac{5}{2}}-2\left (3 a_0+4\right )a_0^{\frac{5}{2}}=0\,,\quad a_w\geqslant a_0\,. \end{equation}

Considering (5.9), the numerical solution for $a_w$ in terms of $a_0$ can be derived, leading to the determination of the critical slope

(5.10) \begin{equation} q_{{cr}}=\frac {\sqrt {a_w-a_0}}{\sqrt {1+a_0}}\sim \sqrt {6\varepsilon }\left (1-\frac {19}{20}\varepsilon +\frac {911}{800}\varepsilon ^2\right )+O(\varepsilon ^3)\,. \end{equation}

Thus, when $\theta _0\gt \arctan q_{{cr}}$ , the soliton dynamics demonstrates the X-shaped pattern rather than the Mach stem, as shown in figure 8, termed the regular reflection. The generation mechanism of the X-shaped pattern is the interaction of two line solitons with an ignorable phase shift. The numerical results for the regular reflection of the reverse bent soliton with $\theta _0=45^\circ$ and $a_0=0.21$ in the initial data (4.1) are depicted in figure 8. The X-shaped soliton pattern is observed to evolve symmetrically, with the two solitons intersecting at the same point as the reverse bent soliton propagates.

Figure 8.

Time evolution of the bent-stem soliton initial data (4.1) with $\theta _0=45^\circ$ and $a_0=0.21$ : (a) $t=0$ , (b) $t=150$ and (c) $t=600$ . The snapshots are shown in a frame of reference moving with $v=-1.1/\cos (45^\circ )\approx -1.5556$ in the $x$ direction.

5.3. Comparison between the BL and KP equations

Through the above results, based on the modulation equations (3.1) for the KP equation, we can conclude similar results for the Mach stem in the reverse bent soliton case. Specifically, the velocity $(s_x,s_y)$ at the discontinuity point $P$ is shown as

(5.11) \begin{equation} \begin{aligned} s_x=\frac {g_0+g_i-g_w}{f_0+f_i-f_w}=\frac {G^{(1)}_0-G^{(1)}_w}{F^{(1)}_0-F^{(1)}_w}=\frac {G^{(1)}_i-G^{(1)}_w}{F^{(1)}_i-F^{(1)}_w}=\frac {G^{(1)}_0-G^{(1)}_i}{F^{(1)}_0-F^{(1)}_i}\,,\\ s_y=\frac {h_0+h_i-h_w}{f_0+f_i-f_w}=\frac {H^{(2)}_0-H^{(2)}_w}{F^{(2)}_0-F^{(2)}_w}=\frac {H^{(2)}_i-H^{(2)}_w}{F^{(2)}_i-F^{(2)}_w}=\frac {H^{(2)}_0-H^{(2)}_i}{F^{(2)}_0-F^{(2)}_i}\,, \end{aligned} \end{equation}

with

(5.12) \begin{equation} \begin{gathered} f=a^{3/2}\,,\quad g=-a^{3/2}+\frac {1}{2}a^{3/2}q^2-\frac {3}{10}a^{5/2}\,,\quad h=-a^{3/2}q\,,\\ F^{(1)}=\frac {1}{q}\,,\quad G^{(1)}=-\frac {a}{2q}-\frac {1}{q}-\frac {q}{2}\,,\\ F^{(2)}=q\,,\quad H^{(2)}=-1-\frac {a}{2}-\frac {q^2}{2}\,. \end{gathered} \end{equation}

The integrability of the KP equation implies the existence of exact Y-shaped soliton solutions, and the amplitude and the slope of each leg are uniquely determined by the phase parameters $\kappa _1$ , $\kappa _2$ and $\kappa _3$ (see Biondini Reference Biondini2007). Without loss of generality, we assume $\kappa _1\lt \kappa _2\lt \kappa _3$ . Thus, we have

(5.13) \begin{equation} \begin{gathered} A_0^{({KP})}=\frac {(\kappa _2-\kappa _1)^2}{6\times 2^{1/3}}\,,\quad A_w^{({KP})}=\frac {(\kappa _3-\kappa _1)^2}{6\times 2^{1/3}}\,,\quad A_i^{({KP})}=\frac {(\kappa _3-\kappa _2)^2}{6\times 2^{1/3}}\,,\\ q_0=-\frac {\kappa _2+\kappa _1}{2^{2/3}}\,,\quad q_w=-\frac {\kappa _3+\kappa _1}{2^{2/3}}\,,\quad q_i=-\frac {\kappa _3+\kappa _2}{2^{2/3}}\,. \end{gathered} \end{equation}

We let $\kappa _3=-\kappa _1$ to ensure $q_w=0$ for the symmetry of the reverse bent soliton case. Replacing $\kappa _1$ and $\kappa _2$ with $A_0^{({KP})}$ and $q_0$ makes it straightforward to infer that

(5.14) \begin{equation} \begin{gathered} A_w^{({KP})}=\left ( \frac {q_0}{\sqrt {3}}+\sqrt {A_0^{({KP})}} \right )^2\,,\quad A_i^{({KP})}=\frac {q_0^2}{3}\,,\quad q_i=-\sqrt {3}\sqrt {A_0^{({KP})}}\,,\\ s_x=-1-\frac {A_0^{({KP})}}{2}=c_w^{({KP})}\,,\quad s_y=\frac {1}{3}\left ( \sqrt {3A_0^{({KP})}}-q_0 \right )\,. \end{gathered} \end{equation}

Thus, $s_y=0$ leads to the critical slope $q_{{cr}}^{({KP})}=\sqrt {3A_0^{({KP})}}$ . The scaling (2.10) converts (5.14) to the data for (1.1), which agree well with the findings presented in Ryskamp et al. (Reference Ryskamp, Hoefer and Biondini2022).

Finally, we can obtain the reduction factor and the critical angle in terms of the initial angle $\theta _0$ and speed $\omega _0$ , i.e.

(5.15) \begin{align} \lambda ^{({KP})}& =\frac {A^{({KP})}_w}{A^{({KP})}_0}=\left (\frac {\sqrt {2}\sin \theta _0}{\sqrt {3}\sqrt {-3-\cos 2\theta _0-4\omega _0\cos \theta _0}}+1\right )^2\,, \nonumber \\ \theta _{{cr}}^{({KP})}& =\arctan q^{({KP})}_{{cr}}=\text{arcsec}\left ( \frac {\sqrt {9\omega _0^2-8}-3\omega _0}{4}\right ). \end{align}

Figure 9.

The same as figure 5 except that $\theta _0\gt 0$ for the Mach reflection of the line solitons.

Figure 9 illustrates the comparison of the wall reduction factor and the Mach stem speed for the KP and BL equations as functions of the initial angle $\theta _0$ and speed $\omega _0$ . This figure deviates from the previous one, denoted as figure 5, in that it specifically considers scenarios where the initial angle $\theta _0\gt 0$ , thus focusing on the Mach reflection of line solitons. In figure 9(a), the wall reduce factor $\lambda$ of the BL equation is consistently greater than that of the KP equation across various initial angles and velocities. The disparity between these coefficients diminishes as the initial angle decreases or the speed approaches unity. Analogously, the variation of the Mach stem speed $c_w$ in figure 9(b) exhibits a similar trend. However, diverging from the observations in figure 5(a), it is discernible from figure 9(a) that the critical angle $\theta _{{cr}}^{({BL})}$ is less than $\theta _{{cr}}^{({KP})}$ , and this discrepancy intensifies with an increase in speed $|\omega _0|$ .

6.1. Subcritical speed

Our concern in this subsection is to investigate the solution behaviour of (2.6) when a localized forcing topography moves with a subcritical speed ( $c_b\lt 1$ ). Thus, when the wave propagates a sufficient distance from the topography, the impact of the topography becomes negligible and can, therefore, be disregarded. Consequently, we are able to conduct a detailed analysis of the Whitham equations (4.2). Here, following the approach of Lee & Grimshaw (Reference Lee and Grimshaw1990), we aim to derive a similarity solution under the assumption that $k=1$ and $a$ is independent of $y$ , leading to

(6.2) \begin{equation} \begin{gathered} q=\pm \frac {y}{\sqrt {f^2-y^2}}\,,\quad a=f^2_t-1\,,\quad \eta =x\pm \sqrt {f^2-y^2}-f+C\,, \end{gathered} \end{equation}

where $C$ is an integration constant, and the function $f(t)$ is subject to the constraint

(6.3) \begin{equation} f(t)=\frac {2a(1+a)^{3/2}(5+3a)}{(12a^2+25a+15)a_t}\,. \end{equation}

Figure 10.

Evolution of the wave generated by the moving topography with $b_0=0.1$ and $c_b=0.8$ at (a–c) $t=50$ and (d–f) $t=450$ . (a,d) Contour plots of the wave, where the black dash–dotted circular arcs with a radius of $t$ and a centre at the origin $(0,0)$ denote the predicted locations of the start-up waves and the black dotted line denotes the linear theory prediction of the wake angle $19.96^\circ$ ; (b,e) $y$ cross-sections of the wave; (c,f) positions of the topography at the corresponding times.

From the above results, one can observe that the wave’s shape is a circle with the centre at $(f(t)-C,0)$ and radius $f(t)$ , while the wave propagates at speed $\partial _tf(t)$ . Furthermore, we assume $a=0$ and $t=0$ , which leads to $f(t)=\pm\, t$ , and a start-up wave initiated by abrupt topography movement is a circle with the centre at $(0,0)$ and radius $t$ . Figure 10 shows the evolution of the wave generated by the moving topography with $b_0=0.1$ and $c_b=0.8$ at $t=50$ and $t=450$ . It is depicted in figure 10 that the movement of the topography generates a series of oscillation waves – the so-called transverse waves, the maximum amplitude of which is of the same order as $b_0$ and gradually diminishes to zero in the $y$ direction. Consequently, the variation of $a$ in the $y$ direction is negligible, which guarantees the validity of dropping the $y$ -dependence in theoretical analyses. As depicted in figures 10(a) and 10(d), the black dash–dotted curves represent the predicted positions of the start-up wave based on the abovementioned discussions. It can be observed that the numerical results are in accordance with the theoretical predictions.

In figure 10(d), it can be seen that the moving topography induces a sequence of trailing waves. By focusing on wave wakes, an approximate analysis can be conducted using the linear theory. Linearizing the BL equation (with a flat bottom) yields

(6.4) \begin{equation} \xi _{tt}-\unicode{x1D6E5} \xi _{tt}-\unicode{x1D6E5} \xi =0\,, \end{equation}

which leads to the linear dispersion relation

(6.5) \begin{equation} -\omega ^2+\boldsymbol{k}^2-\omega ^2\boldsymbol{k}^2=0\,. \end{equation}

Thus, considering the speed $c_b$ of the moving topography, the dispersion relation in the moving frame is derived as follows:

(6.6) \begin{equation} G(k,l)=c_b|\boldsymbol{k}|\cos \phi -\frac {|\boldsymbol{k}|}{\sqrt {1+\boldsymbol{k}^2}}=0\,, \end{equation}

where $\phi$ is the inclination of the normal (direction of $\boldsymbol{k}$ ) to the line of the waves. Hence, the wake angle $\alpha$ denoting the slope of a ray through the topography is (see, Whitham Reference Whitham1999)

(6.7) \begin{equation} \tan \alpha =\frac {\partial _lG(k,l)}{\partial _kG(k,l)}=\frac {c_b^2 \tan \phi }{\left (\tan ^2\phi +1\right )^2-c_b^2}\,. \end{equation}

Clearly, $\alpha \to 0$ as $\phi \to 0$ or $\phi \to \pi /2$ . Hence, the maximum value of $\alpha$ occurs when $\phi \in (0, \pi /2)$ , which can be determined as

(6.8) \begin{equation} \alpha _m=\arctan \frac {3 \sqrt {3} c_b^2 \sqrt {\sqrt {4-3 c_b^2}-1}}{4 \left (2-3c_b^2+\sqrt {4-3 c_b^2}\right )},\quad \text{at}\quad \phi _m=\arctan \sqrt {\frac {\sqrt {4-3c_b^2}-1}{3}}\,. \end{equation}

When $c_b=0$ , it follows that $\alpha _m=0$ . The wake angle gradually expands as $c_b$ increases, reaching $90^\circ$ as $c_b$ approaches unity. The dependence of $\alpha _m$ on $c_b$ is a remarkable wave phenomenon in shallow water, completely different from the deep-water situation, where the famous Kelvin wake angle of $19.47^\circ$ is retained regardless of the speed of the moving object on the free surface. Figure 10(d) presents a comparison between the theoretical predictions and numerical results for the wake angle, depicted by the black dashed lines. It is observable that the trailing waves are all confined within the actualized region, which further confirms the validity of the linear theory.

Figure 11.

Evolution of the wave generated by the moving topography with $b_0=0.1$ and $c_b=1$ at (a–c) $t=50$ and (d–f) $t=1400$ . (a,d) contour plots of the wave, where the black dash–dotted circular arcs with a radius of $t$ and a centre at the origin $(0,0)$ denote the predicted locations of the start-up waves and the black dotted circular arc coarsely fitted as $(x-350)^2+y^2=1100^2$ denotes the locations of the first precursor wave; (b,e) $y$ cross-sections of the wave; (c,f) positions of the topography at the corresponding times.

6.2. Critical speed

In the case of moving topography with critical speed ( $c_b=1$ ), the wave generated by topography will propagate at the same speed as the topography, which is depicted in figure 11. In figures 11(a) and 11(d), it can be observed that the start-up wave is ‘propelled’ a certain distance by the topography due to the consistency in speed between the two, leading to the formation of what is known as the upstream-advancing waves, thereby enabling it to surpass the boundary of the theoretical circle (the black dash–dotted curves). When the start-up wave is counter to the direction of topography movement, it remains unaffected by the influence of the topography. As depicted in figures 11(b) and 11(e), it is evident that the ‘propulsion’ by the topography occurs in the form of waves propagating forward. As time evolves, the distance over which the upstream-advancing waves are propelled increases, and the number of precursor waves also augments. During the propulsion process, the amplitude of the first precursor wave undergoes continuous variation; hence, the impact of $a_t$ in (6.2) cannot be disregarded. The simplification by assuming $a=0$ for approximation does not adequately predict the behaviour. However, by disregarding the local wave characteristics near the topography and instead considering the far-field waveform, it can still be approximated by circular functions (see the black dotted curves in figure 11 d). Considering the wake induced by the topography, shown in figure 11(d), the transverse wave has significantly diminished in intensity, while the divergence wave has gained prominence in contrast to figure 10(d). From (6.8), we can deduce that as $c_b\to 1$ , the wake angle approaches $\pi /2$ , implying that the wake will no longer be confined within a specific angle but will gradually diffuse across the entire plane behind the topography, shown in figure 11(d). This diffusion, however, is a slow process, necessitating extended computational time for its evolution.

Figure 12.

Evolution of the first precursor wave with $c_b=1$ , $b_0=0.1$ and the time interval $\unicode{x1D6E5} t=25$ : (a) variation of its amplitude, where the inset shows the amplitude evolution at the plane $y=0$ over time with the amplitude of the topography $b_0=0.05$ , $0.1$ and $0.15$ ; (b) evolution of its crestlines in the $x{-}y$ plane; (c) evolution of the distance between the crestlines and the origin point $(0,0)$ .

To delve deeper into the dynamics of the first precursor wave, we examine the morphological changes of its crestlines and the amplitude variations along these lines, considering a specific topography amplitude of $b_0=0.1$ . Figure 12(a) illustrates the temporal progression of the first precursor wave’s amplitude, starting from $t=50$ and displaying the results at the time interval of $\unicode{x1D6E5} t=25$ until $t=400$ . This figure includes an inset detailing the amplitude’s evolution at the plane $y=0$ over time for various topography amplitudes, specifically $b_0=0.05$ , $0.1$ and $0.15$ . The amplitude evolution is marked by a rapid increase near the initial time, then reaching the maximum value, after which it undergoes a temperate decay. Concurrently with this decay, energy redistribution occurs, manifesting as a lateral diffusion that leads to a progressive levelling of the crestlines. Comparing the amplitudes of the first precursor waves for the three distinct topography amplitudes, one observes that an increase in the topography amplitude leads to a corresponding increase in the growth rate and peak amplitude of the precursor wave. Additionally, the decay of the wave is accelerated, and the amplitude, once it reaches a steady state, is also notably higher. The spatial evolution of the crestlines of the first precursor wave in the $x{-}y$ plane is shown in figure 12(b), where the crestlines evolve in a manner that resembles a circular arc emanating from the origin $(0,0)$ as time elapses. The evolution of the distance between the crestlines and the origin is depicted in figure 12(c). It can be observed that the distance is almost equal at any given moment, with only a slight protrusion near $y=0$ , and this protrusion becomes increasingly evident as time progresses, which indicates that the peak line is essentially an arc and the main reason for the protrusion is the ‘propulsion’ effect of the topography mentioned above, which causes the centre of the circle to move forward. While the above findings are derived from numerical experiments, the precise patterns of how amplitude and peak lines move with the topography still require further theoretical investigations.

6.3. Supercritical speed

When the topography moves with the supercritical speed ( $c_b\gt 1$ ), the upstream waves disappear; instead, a wave emerges ahead of the topography, moving at a similar speed. This wave, under the ‘propelling’ influence of the topography, gradually narrows, and its amplitude slowly increases until it reaches a steady state, as shown in figures 13(b) and 13(e). In this case, the transverse waves vanish entirely, leaving behind only the distinct strip-like divergence waves, which can be observed in figure 13(d). For the wake angle, the condition $c_b\gt 1$ leads to the emergence of complex values in (6.8), rendering the formula ineffective. Considering the situation from an alternative perspective, the anterior portion of the strip-like divergence wave, influenced by the topography, advances at a speed of $c_b$ , while the rear part of the strip-like divergence wave, which is distant from the topography, is confined to the characteristic wave speed of the fluid, and can only propagate outwards at the unit speed, which can also be discerned from figures 13(a) and 13(d). Therefore, the strip-like divergence wave is confined within a specific angle, which starkly contrasts the wake angle for the topography moving at a subcritical speed. The wake angle $\alpha$ in the supercritical case can be obtained as $\alpha =\arcsin (1/c_b )$ , which (the black dotted line) exhibits a commendable consistency with the numerical simulation shown in figures 13(a) and 13(d).

Figure 13.

Evolution of the wave generated by the moving topography with $b_0=0.1$ and $c_b=1.5$ at (a–c) $t=50$ and (d–f) $t=300$ . (a,d) Contour plots of the wave, where the black dash–dotted circular arcs with a radius of $t$ and a centre at the origin $(0,0)$ denote the predicted locations of the start-up waves and the black dotted lines denote the linear predictions of wake angles $41.81^\circ$ ; (b,e) $y$ cross-sections of the wave; (c,f) the positions of the topography at the corresponding times.

In a realistic situation, the bottom of the ocean may be complicated, with multiple topographies. When a homogeneous or density-stratified flow passes through such a region, the interactions of wave wakes may happen. To investigate these interactions, we use two Gaussian functions to describe the bottom configuration, i.e.

(6.9) \begin{equation} b(x,y,t)=B(x-c_bt,y-l/2)+B(x-c_bt,y+l/2)\,, \end{equation}

where $l$ measures the separation distance between two locally confined topographies in the $y$ direction. The numerical results for the wakes generated by the two topographies moving at the same speed are depicted in figure 14. When $t=400$ , the wakes emanating from the topographies commence to interact with each other. Up to $t=800$ , the wake pattern exhibits a typical reflection characterised by the $y=0$ plane acting as the plane of symmetry. As time passes, the increase of the inner angle between the wakes, which leads to a gradual increase in amplitude, fulfils the conditions conducive to the onset of a Mach reflection. This eventually results in the formation of a Mach stem between the two wakes, as demonstrated in figure 14(c).

Figure 14.

Evolution of the wave generated by the moving topography (6.9) with $b_0=0.1$ , $c_b=1.1$ and $l=400$ at $t=400$ (a), $800$ (b) and $1200$ (c), where the black cross ‘ $\times$ ’ denotes the positions of the topography peaks. The snapshots are shown in a frame of reference moving with $v=1.1$ in the $x$ direction.

We conducted numerical simulations to explore whether altering topography parameters can lead to Mach reflection through wake interactions. The results for various topography speeds and amplitudes are illustrated in figure 15. Our findings indicate that Mach reflection does not occur when the amplitude of the topography is small or when the speed of the topography is large, which corresponds to a smaller angle between the two wakes. In these tests, only regular reflection is observed. Further investigations are needed to determine the specific conditions under which the Mach stem phenomenon occurs, particularly regarding the amplitude and velocity of the topography.

Figure 15.

(a i,b i,c i) Evolution of the wave generated by the moving topography (6.9) with different amplitudes and speeds at $t=1800$ and $l=400$ . (a ii,b ii,c ii) The black solid and dotted curves denote the wave profiles at $y=0$ and $y=200$ , respectively. (a iii,b iii,c iii) The dashed curves denote the positions of the topography at the corresponding parameters.

Note that we fix the distance between two Gaussian topographies in figures 14 and 15, setting $l=400$ . If $l$ is increased, the interaction angle may change, potentially leading to a weak interaction. To investigate the condition under which the Mach stem appears, we conduct numerical experiments by varying $l$ while keeping $b_0=0.1$ and $c_b=1.1$ (the same data as figure 14). In summary, the Mach stem is generated when $l\leqslant 1000$ with the generation time $t\in (0, 3001)$ . Additionally, increasing $l$ results in longer generation times. For instance, we observe the Mach stem for $l=100$ , $200$ and $1000$ when $t\in (0, T]$ , where $T=490$ , $1200$ and $3000.5$ , respectively. However, the situation for $l\gt 1000$ remains to be explored, as it requires significantly more time for numerical simulations, and a thorough analysis will be conducted in the future.