3.1. Leading-order potential

By the double Fourier transforms $\bar {\mathcal {F}}= \int _{-\infty }^{\infty }\mathcal {F}\,\textrm {e}^{\textrm {i}\omega t}\textrm {d} t, \bar {\bar {\mathcal {F}}}= \int _{-\infty }^{\infty }\bar {\mathcal {F}}\,\textrm {e}^{-\textrm {i} kx}\textrm {d}\kern0.06em x$, with $\omega$ the angular velocity and $k$ the wavenumber, equations  (3.5), (3.6) and (3.7) become

(3.13)$$\begin{gather} \frac{\partial^{2}\bar{\bar{\phi}}_{l0}}{\partial z^{2}}+\Bigg(\frac{\omega^{2}}{C_{l}^{2}}-k^{2}\Bigg) \bar{\bar{\phi}}_{l0}=0, \end{gather}$$

(3.14)$$\begin{gather}\frac{\partial^{2} \bar{\bar{\phi}}_{s0}}{\partial z^{2}}+\Bigg(\frac{\omega^{2}}{C_{p}^{2}}-k^{2}\Bigg) \bar{\bar{\phi}}_{s0}=0, \end{gather}$$

(3.15)$$\begin{gather}\frac{\partial^{2}\bar{\bar{\boldsymbol{\varPsi}}}_{0}}{\partial z^{2}}+\left(\frac{\omega^{2}}{C_{s}^{2}}-k^{2}\right) \bar{\bar{\boldsymbol{\varPsi}}}_{0}=0. \end{gather}$$

Let $E_{1}, E_{2}, D_{1}, D_{2}$ be unknowns to be solved for. Then the choice exists either to select $r^{2}=({\omega ^{2}}/{C_{l}^{2}}-k^{2})$, with $r\in \mathbb {R}$, leading to a trial solution for $\bar {\bar {\phi }}_{l0}$ in (3.13) of the form $\bar {\bar {\phi }}_{l0}(z)=E_{1}\cos (rz)+E_{2}\sin (rz)$, or to select $r^{2}=(k^{2}-{\omega ^{2}}/{C_{l}^{2}})$, leading to a trial solution of the form $\bar {\bar {\phi }}_{l0}(z)=E_{1}\cos (\textrm {i} rz)+E_{2}\sin (\textrm {i} rz)$, as in (12c) of Eyov et al. (Reference Eyov, Klar, Kadri and Stiassnie2013). To maintain compatibility with Eyov et al. (Reference Eyov, Klar, Kadri and Stiassnie2013) we choose $r^{2}=(k^{2}-{\omega ^{2}}/{C_{l}^{2}})$. For (3.14) and (3.15) we take $q^{2}=(k^{2}-{\omega ^{2}}/{C_{p}^{2}})$ and $s^{2}=(k^{2}-{\omega ^{2}}/{C_{s}^{2}})$. As in Eyov et al. (Reference Eyov, Klar, Kadri and Stiassnie2013), $r, q$ and $s$ are wavenumbers.

To arrive at a trial solution for (3.14) we choose $\bar {\bar {\phi }}_{s0}(z)=D_{1}\,\textrm {e}^{qz}$, because $\bar {\bar {\phi }}_{s0}(z)\rightarrow 0$ as $z\rightarrow -\infty$ implies no terms involving $\textrm {e}^{-qz}$. By a similar argument $\bar {\bar {\boldsymbol {\varPsi }}}_{0}(z)=D_{2}\,\textrm {e}^{sz}\boldsymbol {j}$.

Applying the double Fourier transforms to the boundary condition at $z=0$ (leading-order term):

(3.16)\begin{equation} \frac{\partial^{2}\phi_{l0}}{\partial t^{2}}+g\frac{\partial\phi_{l0}}{\partial z}=0, \quad z=0. \end{equation}

After both Fourier transforms equation (3.16) becomes

(3.17)\begin{equation} -\omega^{2}\bar{\bar{\phi}}_{l0}+g \frac{\partial\bar{\bar{\phi}}_{l0}}{\partial z}=0, \quad z=0. \end{equation}

Applying Fourier transforms to the first boundary condition at $z=-h$ (2.9), leading-order terms become

(3.18a,b)$$\begin{gather} w_{s0}=\frac{\partial\phi_{s0}}{\partial z}+ \frac{\partial\psi_{0y}}{\partial x} \quad \text{and} \quad \dot{w}_{l0}={-}\frac{\partial\phi_{l0}}{\partial z}, \end{gather}$$

(3.19)$$\begin{gather}-\frac{\partial\phi_{l0}}{\partial z}=\frac{\partial^{2}\phi_{s0}}{\partial t\partial z}+\frac{\partial^{2}\psi_{0y}}{\partial t\partial x}+W(x,y,t). \end{gather}$$

It is only necessary to apply the relevant transforms to the first three terms of (3.19) – the required transforms for $W(x,y,t)$ are already known from Mei & Kadri (Reference Mei and Kadri2018). Assembling terms gives

(3.20)\begin{equation} -\frac{\partial\bar{\bar{\phi}}_{l0}}{\partial z}={-}\textrm{i}\omega\frac{\partial\bar{\bar{\phi}}_{s0}}{\partial z} +\omega k\bar{\bar{\psi}}_{0y}+ \frac{4W_{0}\sin(kb)\sin(\omega T)}{k\omega}, \quad z={-}h. \end{equation}

The second boundary condition at $z=-h$ is $\sigma _{zz}=-P_{l}$:

(3.21)\begin{equation} \sigma_{zz}=\lambda\left(\frac{\partial u_{s0}}{\partial x}+\frac{\partial w_{s0}}{\partial z}\right)+2\mu\frac{\partial w_{s0}}{\partial z}={-}P_{l0}={-}\rho_{l}\frac{\partial\phi_{l0}}{\partial t}, \end{equation}

with

(3.22a,b)\begin{equation} u_{s0}=\frac{\partial\phi_{s0}}{\partial x}-\frac{\partial\psi_{0y}}{\partial z}, \quad \,w_{s0}=\frac{\partial\phi_{s0}}{\partial z}+\frac{\partial\psi_{0y}}{\partial x}. \end{equation}

After application of both Fourier transforms we have

(3.23)\begin{equation} \lambda\left({-}k^{2}\bar{\bar{\phi}}_{s0}+ \frac{\partial^{2}\bar{\bar{\phi}}_{s0}}{\partial z^{2}}\right) +2\mu\left(\frac{\partial^{2}\bar{\bar{\phi}}_{s0}}{\partial z^{2}} +\textrm{i} k\frac{\partial\bar{\bar{\psi}}_{0y}}{\partial z} \right)=\textrm{i}\rho_{l}\omega\bar{\bar{\phi}}_{l0}, \quad z={-}h. \end{equation}

The third boundary condition at $z=-h$ is $\sigma _{xz}=0 \Rightarrow$

(3.24)\begin{equation} \frac{\partial u_{s0}}{\partial z}+\frac{\partial w_{s0}}{\partial x}=0. \end{equation}

After application of both Fourier transforms we have

(3.25)\begin{equation} 2\,\textrm{i} k\frac{\partial\bar{\bar{\phi}}_{s0}}{\partial z}-k^{2}\bar{\bar{\psi}}_{0y}-\frac{\partial^{2} \bar{\bar{\psi}}_{0y}}{\partial z^{2}}=0, \quad z={-}h. \end{equation}

Finally the fourth boundary condition at $z=-h$ is $\sigma _{yz}=0 \Rightarrow$

(3.26)\begin{equation} \frac{\partial v_{s0}}{\partial z}+\frac{\partial w_{s0}}{\partial y}=0, \end{equation}

with

(3.27a,b)\begin{equation} v_{s0}=0 \quad \text{and}\quad w_{s0}=\frac{\partial\phi_{s0}}{\partial z}+\frac{\partial\psi_{0y}}{\partial x}. \end{equation}

Again, after application of both Fourier transforms we arrive at

(3.28)\begin{equation} \frac{\partial^{2}\bar{\bar{\phi}}_{s0}}{\partial z\partial y}+\textrm{i} k\frac{\partial\bar{\bar{\psi}}_{0y}}{\partial y}=0, \quad z={-}h. \end{equation}

3.2. Transformed governing equations

We assemble terms and drop the zero subscript for ease of notation (but remembering these are leading-order terms). Also, we drop the double overbar, again remembering that these are the terms after the double Fourier transforms:

(3.29a–d)$$\begin{gather} \varPhi_{s} = \bar{\bar{\phi}}_{s0}, \quad \varPhi_{l} =\bar{\bar{\phi}}_{l0}, \quad \psi_{y} = \bar{\bar{\psi}}_{0y}, \quad \varPsi = \bar{\bar{\varPsi}}, \end{gather}$$

(3.30a–c)$$\begin{gather}\frac{\partial^{2}\varPhi_{l}}{\partial z^{2}}+\Bigg(\frac{\omega^{2}}{C_{l}^{2}}-k^{2}\Bigg)\varPhi_{l}=0, \enspace \frac{\partial^{2}\varPhi_{s}}{\partial z^{2}}+\Bigg(\frac{\omega^{2}}{C_{p}^{2}}-k^{2}\Bigg)\varPhi_{s}=0, \enspace \frac{\partial^{2}\boldsymbol{\varPsi}}{\partial z^{2}}+\left(\frac{\omega^{2}}{C_{s}^{2}}-k^{2}\right) \boldsymbol{\varPsi}=0. \end{gather}$$

At $z=0$ we have the (transformed) boundary condition for the liquid surface:

(3.31)\begin{equation} -\omega^{2}\varPhi_{l}+g\frac{\partial\varPhi_{l}}{\partial z}=0. \end{equation}

Then, at $z=-h$ we have four (transformed) boundary conditions for the seabed:

(3.32)$$\begin{gather} -\frac{\partial\varPhi_{l}}{\partial z}={-}\textrm{i}\omega\frac{\partial\varPhi_{s}}{\partial z}+\omega k\psi_{y}+\frac{4W_{0}\sin(kb)\sin(\omega T)}{k\omega}, \end{gather}$$

(3.33)$$\begin{gather}\lambda\left({-}k^{2}\varPhi_{s}+\frac{\partial^{2}\varPhi_{s}}{\partial z^{2}}\right)+2\mu\left(\frac{\partial^{2}\varPhi_{s}}{\partial z^{2}}+\textrm{i} k\frac{\partial\psi_{y}}{\partial z}\right)=\textrm{i}\rho_{l}\omega\varPhi_{l}, \end{gather}$$

(3.34)$$\begin{gather}2\,\textrm{i} k\frac{\partial\varPhi_{s}}{\partial z}-k^{2}\psi_{y}-\frac{\partial^{2}\psi_{y}}{\partial z^{2}}=0, \end{gather}$$

(3.35)$$\begin{gather}\frac{\partial^{2}\varPhi_{s}}{\partial z\partial y}+\textrm{i} k\frac{\partial\psi_{y}}{\partial y}=0. \end{gather}$$

With the requirement that $\varPhi _{l}, \varPhi _{s}$ and $\boldsymbol {\varPsi }$, along with all their derivatives, decay to zero as $y\rightarrow \pm \infty, z\rightarrow -\infty$.

3.3. Form for potentials

Taking $\varPhi _{l}(z)=E_{1}\cos (\textrm {i} rz)+E_{2}\sin (\textrm {i} rz)$, we substitute into the boundary condition at $z=0$ to arrive at

(3.36)\begin{equation} E_{2}={-}\frac{\textrm{i}\omega^{2}}{gr}E_{1}, \end{equation}

in agreement with (14a) of Eyov et al. (Reference Eyov, Klar, Kadri and Stiassnie2013). Also, we take $\varPhi _{s}(z)=D_{1}\,\textrm {e}^{qz}+\hat {D}_{1}\,\textrm {e}^{-qz}$ and $\boldsymbol {\varPsi }(z)=\psi _{y}\,\boldsymbol {j}$ with $\psi _{y}=D_{2}\,\textrm {e}^{sz}+\hat {D}_{2}\,\textrm {e}^{-sz}$, but note that, in order to obtain physical solutions in which solid displacements decrease with depth, we must have $\hat {D}_{1}=\hat {D}_{2}=0$ and $s, q \in \mathbb {R}_{\geq 0}$. If this were not the case then displacements would oscillate or increase with depth (Eyov Reference Eyov2013).

Applying boundary condition $\sigma _{xz}=0$ (3.34) at $z=-h$:

(3.37)\begin{equation} D_{2}=\frac{2\,\textrm{i} kq}{k^{2}+s^{2}}\, \textrm{e}^{h(s-q)}D_{1}. \end{equation}

Applying boundary condition (3.32) at $z=-h$:

(3.38)\begin{align} &-E_{1}r\sinh(rh)+\frac{\omega^{2}E_{1}}{g}\cosh(rh)- \textrm{i}\omega qD_{1}\,\textrm{e}^{{-}qh}\nonumber\\ &\quad +\frac{2\,\textrm{i}\omega k^{2}q\, \textrm{e}^{h(s-q)}D_{1}\,\textrm{e}^{{-}sh}}{k^{2}+s^{2}}+ \frac{4W_{0}\sin(kb)\sin(\omega T)}{\omega k}=0. \end{align}

Applying boundary condition (3.33) at $z=-h$:

(3.39)\begin{align} &\lambda({-}k^{2}D_{1}\,\textrm{e}^{{-}qh}+D_{1}q^{2}\,\textrm{e}^{{-}qh}) +2\mu\left(D_{1}q^{2}\,\textrm{e}^{{-}qh}-\frac{2k^{2}D_{1}q\,\textrm{e}^{{-}q h}s}{k^{2}+s^{2}}\right)\nonumber\\ &\quad -\textrm{i}\rho_{l}\omega \left(E_{1}\cosh(rh)-\frac{\omega^{2}E_{1}}{gr}\sinh(rh)\right)=0. \end{align}

Since (3.38) and (3.39) are essentially a pair of simultaneous equations in unknowns $E_{1}$ and $D_{1}$ they can be solved in this case, resulting in

(3.40a,b)\begin{equation} E_{1}={-}\frac{H_{1}}{\omega kH_{2}}, \quad D_{1}=\frac{H_{3}}{kH_{2}}, \end{equation}

with

(3.41)\begin{equation} H_{1}=4grW_{0}\sin(kb)\sin(\omega T)\,\textrm{e}^{{-}qh} \left[{-}2k^{2}\mu qs+(k^{2}+s^{2})\left[\left(\mu+ \frac{\lambda}{2}\right)q^{2}-\frac{\lambda k^{2}}{2}\right]\right] \end{equation}

and

(3.42)\begin{align} H_{2} & =\left({-}2qk^{2}\omega^{2}r\Bigg(\mu s+\frac{\rho_{l}g}{2}\Bigg)+(k^{2}+s^{2}) \omega^{2}r\left(q^{2}\left(\mu+\frac{\lambda}{2}\right)\right.\right.\nonumber\\ &\quad \left.\left.+\,\frac{\rho_{l}gq}{2}- \frac{k^{2}\lambda}{2}\right)\right)\, \textrm{e}^{{-}qh}\cosh(rh) +\left(2qk^{2}\left(\mu gr^{2}s+ \frac{\rho_{l}\omega^{4}}{2}\right)\right.\nonumber\\ &\quad \left.-\,(k^{2}+s^{2})\left(gr^{2}\left(\mu+\frac{\lambda}{2} \right)q^{2}+\frac{\rho_{l}\omega^{4}q}{2}-\frac{gr^{2}k^{2} \lambda}{2}\right)\right)\,\textrm{e}^{{-}qh}\sinh(rh), \end{align}

(3.43)\begin{align} H_{3}&=2\textrm{i}\rho_{l}W_{0}(k^{2}+s^{2}) (\omega^{2}\sinh(rh)-gr\cosh(rh))\sin (kb)\sin(\omega T). \end{align}

Setting $H_2=0$ and rearranging yields

(3.44)\begin{equation} \tanh(rh)=\frac{\dfrac{\omega^{2}}{r}\left\{ \rho_{l}q \dfrac{(k^{2}-s^{2})}{(k^{2}+s^{2})}+\dfrac{1}{g} \left[\dfrac{4k^{2}qs\mu}{(k^{2}+s^{2})}- ((\lambda+2\mu)q^{2}-\lambda k^{2})\right] \right\}}{\dfrac{\omega^{4}q\rho_{l}}{gr^{2}} \dfrac{(k^{2}-s^{2})}{(k^{2}+s^{2})}+ \left[\dfrac{4k^{2}qs\mu}{(k^{2}+s^{2})}-( (\lambda+2\mu)q^{2}-\lambda k^{2})\right]}, \end{equation}

which is the dispersion relation (17) of Eyov et al. (Reference Eyov, Klar, Kadri and Stiassnie2013). The zeros of $H_{2}$ (i.e. dispersion relation solutions) locate the poles for the residue calculations that come later. Therefore, we have

(3.45)$$\begin{gather} \varPhi_{l}(z)={-}\frac{H_{1}}{\omega kH_{2}}\left(\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right), \end{gather}$$

(3.46a,b)$$\begin{gather}\varPhi_{s}(z)= \frac{H_{3}}{kH_{2}}\,\textrm{e}^{qz}, \quad \boldsymbol{\varPsi}=\frac{2\,\textrm{i} q}{k^{2}+s^{2}}\frac{H_{3}}{H_{2}}\, \textrm{e}^{h(s-q)+sz}\boldsymbol{j}. \end{gather}$$

Setting $q=s=0$ reduces to the rigid case where (3.1) of Williams et al. (Reference Williams, Kadri and Abdolali2021) is recovered. Then $H_1$ and $H_2$ reduce to

(3.47)$$\begin{gather} H_{1}={-}2W_{0}gr\lambda k^{4}\sin(kb)\sin(\omega T), \end{gather}$$

(3.48)$$\begin{gather}H_{2}={-}\tfrac{1}{2}r\lambda k^{4}(\omega^{2}\cosh(rh)-gr\sinh(rh)). \end{gather}$$

In this case, $\varPhi _{l}(z)$ becomes

(3.49)\begin{equation} \varPhi_{l}(z)={-}\frac{4W_{0}\sin(kb)\sin(\omega T)}{\mu k\omega}\left\{ \frac{\mu g\cos(\mu z)+\omega^{2}\sin(\mu z)}{\omega^{2}\cos(\mu h)+\mu g\sin(\mu h)}\right\}, \end{equation}

which is in agreement with (3.1) of Williams et al. (Reference Williams, Kadri and Abdolali2021) (note that the sign difference is due to the definition of the velocity potential).

3.4. Inverse Fourier transforms

The leading-order potentials are retrieved by applying the inverse Fourier transforms as

(3.50a–c)\begin{align} \phi_{l0}=\frac{1}{2{\rm \pi}}\int_{-\infty}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}I_{1}, \quad \phi_{s0}=\frac{1}{2{\rm \pi}}\int_{-\infty}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}I_{2}, \quad \psi_{0y}=\frac{1}{2{\rm \pi}}\int_{-\infty}^{\infty}\textrm{i}\,\textrm{d}\omega\,\textrm{e}^{-\textrm{i}\omega t}I_{3}, \end{align}

where $I_{1}, I_{2}, I_{3}$ are the $k$ integrals:

(3.51a,b)$$\begin{gather} I_{1}=\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\,\textrm{e}^{\textrm{i} kx}\frac{-H_{1}}{\omega kH_{2}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right], \quad I_{2}=\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\,\textrm{e}^{\textrm{i} kx}\frac{H_{3}}{kH_{2}}\,\textrm{e}^{qz}, \end{gather}$$

(3.52)$$\begin{gather}I_{3}=\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\, \textrm{e}^{\textrm{i} kx}\frac{2\,\textrm{i} q}{k^{2}+s^{2}}\frac{H_{3}}{H_{2}}\, \textrm{e}^{h(s-q)+sz}. \end{gather}$$

In each case the integrand has poles at the zeros of $H_{2}$, i.e. whenever the dispersion relation (3.44) is satisfied. We substitute out $r, q$ and $s$ to make $H_{2}$ purely a function of $k$. Then values for $I_{1}, I_{2}, I_{3}$ can be calculated from the residues.

Figure 2 shows the various zones where $r, q$ and $s$ take on real and imaginary values. There are zones corresponding to surface waves and acoustic–gravity waves. The remaining zones close to $k=0$ are not physical solutions, since imaginary values taken on by $q$ and/or $s$ would imply oscillations at infinite depth in the elastic medium. Moreover, $q$ and $s$ have to be real and non-negative, otherwise oscillations would increase with increasing depth into the elastic medium. Examination of $I_{1}, I_{2}, I_{3}$ indicates that possible poles might also exist at $k=0$ and when $k^{2}+s^{2}=0$. When $k=0$ the $\sin (kb)$ term in the numerator (from $H_{1}$ and $H_{3}$) ensures a factor of $b$ is reached in the limit $k\rightarrow 0$, so $k=0$ is a removable singularity. For the case $k^{2}+s^{2}=0$ there is a possible pole when $k={\omega }/{\sqrt {2}C_{s}}$, but this pole lies in the unphysical zone of figure 2. From Eyov et al. (Reference Eyov, Klar, Kadri and Stiassnie2013) we have that $s=0$ (at $k_{s}$) represents a point where the energy spreads out over the whole solid depth. At that point the wave amplitude vanishes and so propagation ceases.

Figure 2. Zones possible according to $r, q, s$ being real or imaginary for the case $\omega = 2{\rm \pi}, C_{l}=1450\ {\rm m s}^{-1}, C_{s}=3550\ {\rm m s}^{-1}, C_{p}=6300\ {\rm m s}^{-1}$. Zone 1 (orange) has $r,q,s\in \mathbb {R}$ and corresponds to surface–gravity waves. Zone 2 (green) has $r \in \textrm {i}\mathbb {R}, {\rm with}\ q, s \in \mathbb {R}$ and corresponds to acoustic–gravity waves. The remaining zones near $k=0$ (grey) are not physical solutions. The points where $r,s,q$ transition real $\rightleftharpoons$ imaginary are designated $\pm k_{r}=\pm 0.00433$ (black circles), $\pm k_{s}=\pm 0.00177$ (red circles) and $\pm k_{q}=\pm 0.00099$ (blue circles) respectively.

When $r\Rightarrow r_{0m}$ with $m=0,1$,

(3.53a,b)\begin{equation} r=\sqrt{k^{2}-{\omega^{2}}/{C_{l}^{2}}}, \quad \implies k_{0m}=\sqrt{\frac{\omega^{2}}{C_{l}^{2}}+r_{0m}^{2}}, \end{equation}

which corresponds to surface waves. There are two possible modes for surface waves. Mode 00 can propagate if $\omega > \omega _{00}$ – the cut-off frequency for this mode. Mode 01 is the usual tsunami. If instead $r\Rightarrow \textrm {i} r_{n}$, then acoustic–gravity waves are possible, and

(3.54)\begin{equation} k_{n}=\sqrt{\frac{\omega^{2}}{C_{l}^{2}}-r_{n}^{2}}, \end{equation}

up to a maximum value of $n=N$, after which the evanescent waves exist with wavenumber $\varLambda _{n}$:

(3.55)\begin{equation} k_{n}=\textrm{i}\varLambda_{n}=\textrm{i}\sqrt{r_{n}^{2}- \frac{\omega^{2}}{C_{l}^{2}}}=\sqrt{\frac{\omega^{2}}{C_{l}^{2}}-r_{n}^{2}}. \end{equation}

Solutions to the dispersion relation involving acoustic–gravity waves for the case $\omega =2{\rm \pi}$ occur between $k_{s}=0.00177$ and $k_{r}=0.00433$. They are marked with blue diamonds in figure 3. Considering the liquid terms first, we break $\phi _{l0}$ into the different regions according to varying $\omega$. For $r\in \textrm {i}\mathbb {R}$:

(3.56)\begin{align} \phi_{l0}&=\frac{1}{2{\rm \pi}}\int_{-\infty}^{-\omega_{n}}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\,\textrm{e}^{\textrm{i} kx}\frac{-H_{1}}{\omega kH_{2}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right] \nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{-\omega_{n}}^{\omega_{n}}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\,\textrm{e}^{\textrm{i} kx}\frac{-H_{1}}{\omega kH_{2}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{\omega_{n}}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\,\textrm{e}^{\textrm{i} kx}\frac{-H_{1}}{\omega kH_{2}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right], \end{align}

whereas for $r\in \mathbb {R}$:

(3.57)\begin{align} \phi_{l0}&=\frac{1}{2{\rm \pi}}\int_{-\infty}^{0}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\,\textrm{e}^{\textrm{i} kx}\frac{-H_{1}}{\omega kH_{2}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{0}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\,\textrm{e}^{\textrm{i} kx}\frac{-H_{1}}{\omega kH_{2}}\left[\cos(\textrm{i} rz)- \frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{-\infty}^{-\omega_{00}}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\,\textrm{e}^{\textrm{i} kx}\frac{-H_{1}}{\omega kH_{2}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{\omega_{00}}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\,\textrm{e}^{\textrm{i} kx}\frac{-H_{1}}{\omega kH_{2}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]. \end{align}

Application of the Rayleigh damping method and contour integration using the residue theorem around a positively oriented simple closed curve as per Mei & Kadri (Reference Mei and Kadri2018) results in

(3.58)\begin{align} \phi_{l0}&=\frac{1}{2{\rm \pi}}\int_{-\infty}^{-\omega_{n}}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\sum_{n=1}^{N}\frac{-H_{1}|_{k={-}k_{n}}}{-\omega k_{n}\partial_{k}H_{2}|_{k={-}k_{n}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{-\textrm{i} k_{n}x}\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{-\omega_{n}}^{\omega_{n}}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\sum_{n=N+1}^{\infty}\frac{-H_{1}|_{k=\textrm{i}\varLambda_{n}}}{\omega \textrm{i}\varLambda_{n}\partial_{k}H_{2}|_{k=\textrm{i}\varLambda_{n}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{-\varLambda_{n}x}\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{\omega_{n}}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\sum_{n=1}^{N}\frac{-H_{1}|_{k=k_{n}}}{\omega k_{n}\partial_{k}H_{2}|_{k=k_{n}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{\textrm{i} k_{n}x}, \end{align}

with $r=\textrm {i} r_{n}\in \textrm {i}\mathbb {R}$, $\partial _{k}={\partial }/{\partial k}$,

(3.59)\begin{gather} H_{1}|_{k=k_{n}}=4\textrm{i} gr_{n}W_{0}\sin(k_{n}b)\sin(\omega T)\,\textrm{e}^{{-}q_{n}h}\left[{-}2k_{n}^{2}\mu q_{n}s_{n}\vphantom{\frac{\lambda}{2}}\right.\nonumber\\ \quad\qquad \ \ \left.+\,(k_{n}^{2}+s_{n}^{2})\left[\left(\mu+ \frac{\lambda}{2}\right)q_{n}^{2}-\frac{\lambda k_{n}^{2}}{2}\right]\right], \end{gather}

(3.60)\begin{gather} H_{1}|_{k={-}k_{n}}={-}H_{1}|_{k=k_{n}}, \end{gather}

(3.61)\begin{gather} H_{1}|_{k=\textrm{i}\varLambda_{n}}=4\textrm{i} g r_{n}W_{0}\sin(\textrm{i}\varLambda_{n}b)\sin(\omega T)\,\textrm{e}^{{-}q_{n}h}\left[2\varLambda_{n}^{2}\mu q_{n}s_{n}\vphantom{\frac{\lambda}{2}}\right.\nonumber\\ \qquad\qquad\ \left.+\,(s_{n}^{2}-\varLambda_{n}^{2}) \left[\left(\mu+\frac{\lambda}{2}\right)q_{n}^{2}+ \frac{\lambda\varLambda_{n}^{2}}{2}\right]\right]. \end{gather}

The derivative terms are given in Appendix A.

Figure 3. Acoustic–gravity wave solutions to the dispersion relation are located at the intersections of dashed and solid curves (blue diamonds) for $\omega =2{\rm \pi}$ and depth $h=4000$ m. Dashed curve is left-hand side (LHS) of (3.44) and solid curve is right-hand side (RHS) of (3.44) when $r\in \textrm {i}\mathbb {R}$.

In support of the validity of the integration process figure 4 shows a plot of ${1}/{|H_{2}|}$ in the complex plane when $H_{2}=H_{2}(k)$ and $k$ is allowed to take on complex values. Cross-sections through the real and imaginary axes appear in figures 5(a) and 5(b) respectively. The poles of the function lie on the real axis, whereas the zeros lie on the imaginary axis. If the range of the plots were to be extended then the function decays to zero everywhere. As empirical evidence for the validity of the integration, when the calculations are complete, we find good agreement with existing synthetic and real data plots for both acoustic–gravity waves and surface waves (e.g. see figures 13 and 17).

Figure 4. Plot of ${1}/{|H_{2}|}$ in the complex plane when $H_{2}=H_{2}(k)$ and $k$ is allowed to take on complex values. The angular frequency in this case is $\omega =2{\rm \pi}$ as in figure 3.

Figure 5. (a) Cross-section of figure 4 through the real axis showing locations of the poles when $\omega =2{\rm \pi}$. (b) Cross-section of figure 4 through the imaginary axis showing locations of the zeros when $\omega =2{\rm \pi}$.

In the case where $r\Rightarrow r_{0m}, k\Rightarrow k_{0m}$ with $r_{0m}$ a real number and $m=0,1$ then there may exist two possibilities for surface waves:

(3.62) \begin{align} \phi_{l0}&=\frac{1}{2{\rm \pi}}\int_{-\infty}^{0}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{-H_{1}|_{k={-}k_{01}}}{-\omega k_{01}\partial_{k}H_{2}|_{k={-}k_{01}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr_{01}}\sin(\textrm{i} rz)\right]\textrm{e}^{-\textrm{i} k_{01}x} \nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{0}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{-H_{1}|_{k=k_{01}}}{\omega k_{01}\partial_{k}H_{2}|_{k=k_{01}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr_{01}}\sin(\textrm{i} rz)\right]\textrm{e}^{\textrm{i} k_{01}x}\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{-\infty}^{-\omega_{00}}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{-H_{1}|_{k={-}k_{00}}}{-\omega k_{00}\partial_{k}H_{2}|_{k={-}k_{00}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr_{00}}\sin(\textrm{i} rz)\right]\textrm{e}^{-\textrm{i} k_{00}x}\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{\omega_{00}}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{-H_{1}|_{k=k_{00}}}{\omega k_{00}\partial_{k}H_{2}|_{k=k_{00}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr_{00}}\sin(\textrm{i} rz)\right]\textrm{e}^{\textrm{i} k_{00}x}, \end{align}

with

(3.63)\begin{align} H_{1}|_{k=k_{0m}}&=4gr_{0m}W_{0}\sin(k_{0m}b)\sin(\omega T)\,\textrm{e}^{{-}q_{0m}h}\left(\vphantom{\left.+(k_{0m}^{2}+s_{0m}^{2})\left(q_{0m}^{2} \left(\mu+\frac{\lambda}{2}\right)-\frac{\lambda k_{0m}^{2}}{2}\right)\right)}{-}2k_{0m}^{2}\mu q_{0m}s_{0m}\right.\nonumber\\ &\quad \left.+\,(k_{0m}^{2}+s_{0m}^{2})\left(q_{0m}^{2} \left(\mu+\frac{\lambda}{2}\right)-\frac{\lambda k_{0m}^{2}}{2}\right)\right), \end{align}

(3.64)\begin{equation} H_{1}|_{k={-}k_{0m}}={-}H_{1}|_{k=k_{0m}}. \end{equation}

The derivative term is again to be found in Appendix A. Using the substitutions

(3.65a–c)\begin{equation} k=\sqrt{r^2+\frac{\omega^2}{C_{l}^2}},\quad q=\sqrt{r^2+\frac{\omega^2}{C_{l}^2} - \frac{\omega^2}{C_{p}^2}},\quad s=\sqrt{r^2+\frac{\omega^2}{C_{l}^2} - \frac{\omega^2}{C_{s}^2}}, \end{equation}

the dispersion relation (3.44) can be written in terms of $r$ and $\omega$ alone, and in this case the condition for the existence of the $00$th mode for a particular $\omega$ is

(3.66) \begin{equation} \frac{\textrm{d}}{\textrm{d} r} \bigl[{\rm left-hand side of } (3.44)\bigr] > \frac{\textrm{d}}{\textrm{d} r} \bigl[{\rm right-hand side of } (3.44)\bigr]. \end{equation}

The expressions for the velocity potential in the liquid layer can be further reduced to

(3.67)\begin{align} \phi_{l0}&=\frac{1}{\rm \pi}\int_{\omega_{n}}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\sum_{n=1}^{N}\frac{-H_{1}|_{k=k_{n}}}{\omega k_{n}\partial_{k}H_{2}|_{k=k_{n}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{\textrm{i} k_{n}x}\nonumber\\ &\quad +\frac{1}{\rm \pi}\int_{0}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{-H_{1}|_{k=k_{01}}}{\omega k_{01}\partial_{k}H_{2}|_{k=k_{01}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{\textrm{i} k_{01}x}\nonumber\\ &\quad +\frac{1}{\rm \pi}\int_{\omega_{00}}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{-H_{1}|_{k=k_{00}}}{\omega k_{00}\partial_{k}H_{2}|_{k=k_{00}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{\textrm{i} k_{00}x}\nonumber\\ &\quad +\frac{1}{\rm \pi}\int_{0}^{\omega_{n}}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\sum_{n=n+1}^{\infty}\frac{-H_{1}|_{k=\textrm{i}\varLambda_{n}}}{\omega \textrm{i}\varLambda_{n}\partial_{k}H_{2}|_{k=\textrm{i}\varLambda_{n}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{-\varLambda_{n}x}. \end{align}

Returning to the expression for the displacement potential in the solid given by

(3.68)\begin{equation} \phi_{s0}=\frac{1}{2{\rm \pi}}\int_{-\infty}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\,\textrm{e}^{\textrm{i} kx}\frac{H_{3}}{kH_{2}}\,\textrm{e}^{qz}, \end{equation}

and following a similar procedure to that of the liquid velocity potential case, we arrive at

(3.69)\begin{align} \phi_{s0}& =\frac{1}{\rm \pi}\int_{\omega_{n}}^{\infty}\textrm{i}\,\textrm{d}\omega\textrm{e}^{-\textrm{i}\omega t}\sum_{n=1}^{N}\frac{H_{3}|_{k=k_{n}}}{k_{n}\partial_{k}H_{2}|_{k=k_{n}}}\, \textrm{e}^{q_{n}z}\, \textrm{e}^{\textrm{i} k_{n}x}\nonumber\\ &\quad +\frac{1}{\rm \pi}\int_{0}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{H_{3}|_{k=k_{01}}}{k_{01}\partial_{k}H_{2}|_{k=k_{01}}}\, \textrm{e}^{q_{01}z}\, \textrm{e}^{\textrm{i} k_{01}x}\nonumber\\ &\quad +\frac{1}{\rm \pi}\int_{\omega_{00}}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{H_{3}|_{k=k_{00}}}{k_{00}\partial_{k}H_{2}|_{k=k_{00}}}\, \textrm{e}^{q_{00}z}\, \textrm{e}^{\textrm{i} k_{00}x}\nonumber\\ &\quad +\frac{1}{\rm \pi}\int_{0}^{\omega_{n}}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\sum_{n=N+1}^{\infty}\frac{H_{3}|_{k=\textrm{i}\varLambda_{n}}}{\textrm{i}\varLambda_{n}\partial_{k}H_{2}|_{k=\textrm{i}\varLambda_{n}}}\, \textrm{e}^{q_{n}z-\varLambda_{n}x}, \end{align}

with

(3.70) \begin{equation} \left.\begin{array}{c@{}} H_{3}|_{k=k_{n}}={-}2 \rho_{l}W_{0}(k_{n}^{2}+s_{n}^{2}) \sin(k_{n}b)\sin(\omega T)(\omega^{2}\sin(r_{n}h)-gr_{n}\cos(r_{n}h)),\\ H_{3}|_{k=k_{0m}}=2i\rho_{l}W_{0}(k_{0m}^{2}+s_{0m}^{2}) \sin(k_{0m}b)\sin(\omega T)(\omega^{2}\sinh(r_{0m}h)-gr_{0m}\cosh(r_{0m}h)),\\ H_{3}|_{k=\textrm{i}\varLambda_{n}}=2\textrm{i}\rho_{l}W_{0} (s_{n}^{2}-\varLambda_{n}^{2})\sinh (\varLambda_{n}b)\sin(\omega T)(\omega^{2}\sin(r_{n}h)-gr_{n}\cos(r_{n}h)). \end{array}\right\} \end{equation}

The derivative terms evaluated at $k= k_{n}$, $k=\textrm {i}\varLambda _{n}$ and $k=k_{0m}$ remain as before. In a similar way the rotation potential can be written as

(3.71)\begin{equation} \psi_{0y}=\frac{1}{2{\rm \pi}}\int_{-\infty}^{\infty}\textrm{i}\,\textrm{d}\omega\,\textrm{e}^{-\textrm{i}\omega t}\frac{1}{2{\rm \pi} \textrm{i}}\int_{-\infty}^{\infty}\textrm{d} k\,\textrm{e}^{\textrm{i} kx}\frac{2\,\textrm{i} q}{k^{2}+s^{2}}\frac{H_{3}}{H_{2}}\,\textrm{e}^{h(s-q)+sz}, \end{equation}

which becomes

(3.72)\begin{align} \psi_{0y}& =\frac{1}{\rm \pi}\int_{\omega_{n}}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\sum_{n=1}^{N}\frac{2\,\textrm{i} q_{n}}{k_{n}^{2}+s_{n}^{2}}\frac{H_{3}|_{k=k_{n}}}{\partial_{k}H_{2}|_{k=k_{n}}}\, \textrm{e}^{h(s_{n}-q_{n})+s_{n}z}\, \textrm{e}^{\textrm{i} k_{n}x}\nonumber\\ & \quad +\frac{1}{\rm \pi}\int_{0}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{2\,\textrm{i} q_{01}}{k_{01}^{2}+s_{01}^{2}}\frac{H_{3}|_{k=k_{01}}}{\partial_{k}H_{2}|_{k=k_{01}}}\, \textrm{e}^{h(s_{01}-q_{01})+s_{01}z}\, \textrm{e}^{\textrm{i} k_{01}x}\nonumber\\ & \quad +\frac{1}{\rm \pi}\int_{\omega_{00}}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\frac{2\,\textrm{i} q_{00}}{k_{00}^{2}+s_{00}^{2}}\frac{H_{3}|_{k=k_{00}}}{\partial_{k}H_{2}|_{k=k_{00}}}\, \textrm{e}^{h(s_{00}-q_{00})+s_{00}z}\, \textrm{e}^{\textrm{i} k_{00}x}\nonumber\\ & \quad +\frac{1}{\rm \pi}\int_{0}^{\omega_{n}}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\sum_{n=N+1}^{\infty}\frac{2\,\textrm{i} q_{n}}{s_{n}^{2}-\varLambda_{n}^{2}}\frac{H_{3}|_{k=\textrm{i}\varLambda_{n}}}{\partial_{k}H_{2}|_{k=\textrm{i}\varLambda_{n}}}\, \textrm{e}^{h(s_{n}-q_{n})+s_{n}z-\varLambda_{n}x}. \end{align}

3.5. Long-range modulation: liquid layer

We introduce unknown envelope factors for the liquid layer $A_{n}^{\pm }(X,Y)$ and $A_{0m}^{\pm }(X,Y)$:

(3.73)\begin{align} \phi_{l0}&=\frac{1}{2{\rm \pi}}\int_{\omega_{n}}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\sum_{n=1}^{N}A_{n}^+\frac{-H_{1}|_{k=k_{n}}}{\omega k_{n}\partial_{k}H_{2}|_{k=k_{n}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr_{n}}\sin(\textrm{i} rz)\right]\textrm{e}^{\textrm{i} k_{n}x}\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{0}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}A_{01}^+\frac{-H_{1}|_{k=k_{01}}}{\omega k_{01}\partial_{k}H_{2}|_{k=k_{01}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{\textrm{i} k_{01}x}\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{\omega_{00}}^{\infty}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}A_{00}^+\frac{-H_{1}|_{k=k_{00}}}{\omega k_{00}\partial_{k}H_{2}|_{k=k_{00}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{\textrm{i} k_{00}x}\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{-\infty}^{-\omega_{n}}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\sum_{n=1}^{N}A_{n}^{-}\frac{-H_{1}|_{k={-}k_{n}}}{-\omega k_{n}\partial_{k}H_{2}|_{k={-}k_{n}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{-\textrm{i} k_{n}x}\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{-\infty}^{0}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}A_{01}^{-}\frac{-H_{1}|_{k={-}k_{01}}}{-\omega k_{01}\partial_{k}H_{2}|_{k={-}k_{01}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{-\textrm{i} k_{01}x}\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\int_{-\infty}^{-\omega_{00}}\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}A_{00}^{-}\frac{-H_{1}|_{k={-}k_{00}}}{-\omega k_{00}\partial_{k}H_{2}|_{k={-}k_{00}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{-\textrm{i} k_{00}x}\nonumber\\ &\quad +\frac{1}{2{\rm \pi}}\left[\int_{-\omega_{n}}^{0}+\int_{0}^{\omega_{n}}\right]\textrm{i}\,\textrm{d}\omega \, \textrm{e}^{-\textrm{i}\omega t}\sum_{n=N+1}^{\infty}\frac{-H_{1}|_{k=\textrm{i}\varLambda_{n}}}{\omega \textrm{i}\varLambda_{n}\partial_{k}H_{2}|_{k=\textrm{i}\varLambda_{n}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{-\varLambda_{n}x}. \end{align}

The initial conditions are given by

(3.74a,b)\begin{equation} A_{n}^{{\pm}}=\begin{cases} 1 & |Y|< {\unicode{x00A3}}=\epsilon L\\ 0 & |Y|> {\unicode{x00A3}}=\epsilon L , \end{cases} \quad A_{0m}^{{\pm}}=\begin{cases} 1 & Y|< {\unicode{x00A3}}=\epsilon L\\ 0 & Y|> {\unicode{x00A3}}=\epsilon L , \end{cases} \quad X=\epsilon^{2}x\rightarrow 0. \end{equation}

The waves are required to vanish far away from, and be symmetric about, the central axis (Mei & Kadri Reference Mei and Kadri2018):

(3.75)\begin{equation} A_{n}^{{\pm}}=A_{0m}^{{\pm}}=0,\quad |Y|\rightarrow\infty,\quad \frac{\partial A_{n}^{{\pm}}}{\partial Y}=\frac{\partial A_{0m}^{{\pm}}}{\partial Y}=0,\quad Y=0. \end{equation}

Proceeding with acoustic modes, taking the time Fourier transform of (3.8) and separating $\bar {\phi }_{l2}$ into three ranges yields

(3.76) \begin{equation} \bar{\phi}_{l2}=\begin{cases} \bar{\phi}_{l2}^+ & \omega_{n}<\omega<\infty\\ \bar{\phi}_{l2}^{e} & -\omega_{n}<\omega<\omega_{n}\\ \bar{\phi}_{l2}^{-} & -\infty<\omega<{-}\omega_{n} . \end{cases} \end{equation}

In the range $\omega _{n}<\omega <\infty$,

(3.77)\begin{equation} \frac{\partial^{2}\bar{\phi}_{l2}^+}{\partial x^{2}}+\frac{\partial^{2}\bar{\phi}_{l2}^+}{\partial z^{2}}+\frac{\omega^{2}}{C_{l}^{2}}\bar{\phi}_{l2}^+{=}-\frac{\partial^{2}\bar{\phi}_{l0}}{\partial Y^{2}}-2\frac{\partial^{2}\bar{\phi}_{l0}}{\partial x\partial X}. \end{equation}

From this point the solution proceeds in an analogous way to that derived by Mei & Kadri (Reference Mei and Kadri2018):

(3.78)\begin{align} &\frac{\partial^{2}\bar{\phi}_{l2}^+}{\partial x^{2}}+\frac{\partial^{2}\bar{\phi}_{l2}^+}{\partial z^{2}}+\frac{\omega^{2}}{C_{l}^{2}}\bar{\phi}_{l2}^+{=}-\textrm{i}\sum_{n=1}^{N}\left[\frac{\partial^{2}A_{n}^+}{\partial Y^{2}}+2\,\textrm{i} k_{n}\frac{\partial A_{n}^+}{\partial X}\right]\nonumber\\ &\quad \frac{-H_{1}|_{k=k_{n}}}{\omega k_{n}\partial_{k}H_{2}|_{k=k_{n}}}\left[\cos(\textrm{i} rz)-\frac{\textrm{i}\omega^{2}}{gr}\sin(\textrm{i} rz)\right]\textrm{e}^{\textrm{i} k_{n}x}. \end{align}

Assuming $\bar {\phi }_{l2}^+$ has solutions in the form $\sum _{n=1}^{N}\xi _{n}^+(\omega,z)\,\textrm {e}^{\textrm {i} k_{n}x}$, then substituting into (3.78) gives

(3.79)\begin{align} &-\sum_{n=1}^{N}\xi_{n}^+k_{n}^{2}\, \textrm{e}^{\textrm{i} k_{n}x}+\sum_{n=1}^{N}\frac{\partial^{2}\xi_{n}^+}{\partial z^{2}}\, \textrm{e}^{\textrm{i} k_{n}x}+\frac{\omega^{2}}{C_{l}^{2}}\sum_{n=1}^{N}\xi_{n}^+\, \textrm{e}^{\textrm{i} k_{n}x}\nonumber\\ &\quad ={-}\textrm{i}\sum_{n=1}^{N}\left[\frac{\partial^{2}A_{n}^+}{\partial Y^{2}}+2\,\textrm{i} k_{n}\frac{\partial A_{n}^+}{\partial X}\right]\frac{-H_{1}|_{k=k_{n}}}{\omega k_{n}\partial_{k}H_{2}|_{k=k_{n}}}F_{n}(z)\,\textrm{e}^{\textrm{i} k_{n}x}. \end{align}

Equating coefficients gives

(3.80)\begin{equation} \frac{\partial^{2}\xi_{n}^+}{\partial z^{2}}+\left(\frac{\omega^2}{C_{l}^2}-k_{n}^2\right) \xi_{n}^+{=}-\textrm{i}\left[\frac{\partial^{2}A_{n}^+}{\partial Y^{2}}+2\,\textrm{i} k_{n}\frac{\partial A_{n}^+}{\partial X}\right]\frac{-H_{1}|_{k=k_{n}}}{\omega k_{n}\partial_{k}H_{2}|_{k=k_{n}}}F_{n}(z), \end{equation}

where

(3.81) \begin{equation} r^2=k^2-\frac{\omega^2}{C_{l}^2} \begin{cases} r=r_{0m},\quad \dfrac{\omega^2}{C_{l}^2}-k_{0m}^2={-}r_{0m}^2,\quad\textrm{surface waves},\\ r=\textrm{i} r_{n},\quad \dfrac{\omega^2}{C_{l}^2}-k_{n}^2={+}r_{n}^2,\quad\textrm{acoustic--gravity waves}, \end{cases} \end{equation}

resulting in

(3.82)\begin{equation} \frac{\partial^{2}\xi_{n}^+}{\partial z^{2}}+r_{n}^{2}\xi_{n}^+{=}-\textrm{i}\left[\frac{\partial^{2}A_{n}^+}{\partial Y^{2}}+2\,\textrm{i} k_{n}\frac{\partial A_{n}^+}{\partial X}\right]\frac{-H_{1}|_{k=k_{n}}}{\omega k_{n}\partial_{k}H_{2}|_{k=k_{n}}}F_{n}(z), \end{equation}

where

(3.83)\begin{equation} F_{n}(z)=\cos({r_{n}z})+\frac{\omega^2}{g r_{n}}\sin({r_{n}z}). \end{equation}

The ground motion is captured at $O(\epsilon ^0$), so the boundary conditions on $\bar {\phi }_{l2}^+$ (and therefore $\xi _{n}^+$) at O($\epsilon ^2$) are

(3.84)$$\begin{gather} -\omega^2\bar{\phi}_{l2}^+{+}g\frac{\partial \bar{\phi}_{l2}^+}{\partial z}=0,\quad z=0, \end{gather}$$

(3.85)$$\begin{gather}\frac{\partial\bar{\phi}_{l2}^+}{\partial z}=0,\quad z={-}h. \end{gather}$$

Here $F_{n}(z)$ is a solution to the boundary value problem:

(3.86)$$\begin{gather} \frac{\partial^{2} F_{n}}{\partial z^{2}}+r_{n}^{2} F_{n} = 0, \end{gather}$$

(3.87)$$\begin{gather}F_{n} = 1,\quad \frac{\partial F_{n}}{\partial z}=\frac{\omega^2}{g},\quad z=0, \end{gather}$$

(3.88)$$\begin{gather}F_{n}=\cos({r_{n}h})-\frac{\omega^2}{g r_{n}}\sin({r_{n}h}),\quad\frac{\partial F_{n}}{\partial z} = r_{n} \sin(r_{n} h ) + \frac{\omega^2}{g}\cos({ r_{n} h}), \quad z={-}h. \end{gather}$$

A similar process could be carried out for surface waves.

The next step is to extract the Schrödinger equation from (3.82) by applying the following Green's identity over the range $-h \leq z \leq 0$:

(3.89)\begin{align} \int_{{-}h}^{0}\left[F_{n}\left(\frac{\partial^{2}\xi_{n}^+}{\partial z^{2}}+r_{n}^{2}\xi_{n}^+\right)-\xi_{n}^+\left(\frac{\partial^{2}F_{n}}{\partial z^{2}}+r_{n}^{2}F_{n}\right)\right]\textrm{d} z=\left[F_{n}\frac{\partial\xi_{n}^+}{\partial z}-\xi_{n}^+\frac{\partial F_{n}}{\partial z}\right]_{{-}h}^{0}. \end{align}

As in Mei & Kadri (Reference Mei and Kadri2018), the result is the Schrödinger equation for the two-dimensional evolution of the envelope factors:

(3.90)\begin{equation} \frac{\partial^{2}A_{n}^+}{\partial Y^{2}}+2\,\textrm{i} k_{n}\frac{\partial A_{n}^+}{\partial X}=0,\quad r\in \textrm{i}\mathbb{R}. \end{equation}

Having obtained the Schrödinger equation (3.90) for the acoustic–gravity wave case in the liquid layer the solution is analogous to that found in Mei & Kadri (Reference Mei and Kadri2018), but with mode properties now incorporating elasticity, via $k_{n}$. The envelope solution is therefore stated here as

(3.91)\begin{align} A_{n}&=\frac{1-{\rm i}}{2}\left\{ C\left(\sqrt{\frac{2}{{\rm \pi}\chi_{n}}}\mathcal{Y_+}\right)+ C\left(\sqrt{\frac{2}{{\rm \pi}\chi_{n}}}\mathcal{Y_{-}}\right)\right\}\nonumber\\ &\quad +\frac{1+{\rm i}}{2}\left\{ S\left(\sqrt{\frac{2}{{\rm \pi}\chi_{n}}}\mathcal{Y_+}\right)+S\left(\sqrt{\frac{2}{{\rm \pi}\chi_{n}}}\mathcal{Y_{-}}\right)\right\},\quad \chi_{n}=\frac{X}{2k_{n}},\quad \mathcal{Y}_{{\pm}}=\frac{{\unicode{x00A3}}\pm Y}{2}, \end{align}

where $C(z)$ and $S(z)$ are Fresnel integrals. A similar process beginning at (3.76) can be applied to derive the expressions for the surface wave mode 01 and mode 00. Finally, the pressure obtained from the velocity potential (3.67) in the liquid (propagating parts) inclusive of envelope factors is given by

(3.92)\begin{align} P&=\frac{\rho_{l}}{\rm \pi}\int_{\omega_{n}}^{\infty}\,\textrm{d}\omega \sum_{n=1}^{N}\frac{-H_{1}|_{k=k_{n}} A_{n}}{k_{n}\partial_{k}H_{2}|_{k=k_{n}}}\left[\cos(rz) +\frac{\omega^{2}}{gr}\sin(rz)\right]\textrm{e}^{\textrm{i} (k_{n}x-\omega t)}\nonumber\\ &\quad +\frac{\rho_{l}}{\rm \pi}\int_{0}^{\infty}\,\textrm{d}\omega \frac{-H_{1}|_{k=k_{01}} A_{01}}{k_{01}\partial_{k}H_{2}|_{k=k_{01}}} \left[\cosh(rz)+\frac{\omega^{2}}{gr}\sinh( rz)\right]\textrm{e}^{\textrm{i}(k_{01}x-\omega t)}\nonumber\\ &\quad +\frac{\rho_{l}}{\rm \pi}\int_{\omega_{00}}^{\infty}\,\textrm{d}\omega \frac{-H_{1}|_{k=k_{00}} A_{00}}{k_{00}\partial_{k}H_{2}|_{k=k_{00}}} \left[\cosh(rz)+\frac{\omega^{2}}{gr}\sinh( rz)\right]\textrm{e}^{\textrm{i}(k_{00}x-\omega t)}. \end{align}

Similarly the surface elevation is given by

(3.93)\begin{align} \eta&=\frac{1}{g{\rm \pi}}\int_{\omega_{n}}^{\infty}\,\textrm{d}\omega \sum_{n=1}^{N}\frac{-H_{1}|_{k=k_{n}} A_{n}}{k_{n}\partial_{k}H_{2}|_{k=k_{n}}}\, \textrm{e}^{\textrm{i} (k_{n}x-\omega t)} +\frac{1}{g{\rm \pi}}\int_{0}^{\infty}\,\textrm{d}\omega \frac{-H_{1}|_{k=k_{01}} A_{01}}{k_{01}\partial_{k}H_{2}|_{k=k_{01}}}\, \textrm{e}^{\textrm{i}(k_{01}x-\omega t)}\nonumber\\ &\quad +\frac{1}{g{\rm \pi}}\int_{\omega_{00}}^{\infty}\,\textrm{d}\omega \frac{-H_{1}|_{k=k_{00}} A_{00}}{k_{00}\partial_{k}H_{2}|_{k=k_{00}}}\, \textrm{e}^{\textrm{i}(k_{00}x-\omega t)}. \end{align}

4.1. Acoustic–gravity waves

When the acoustic–gravity wave propagating modes ($n=2,3,\ldots$) terminate, the phase velocity becomes equal to $C_{s}$ and $s=\sqrt {k^2-{\omega ^2}/{C_{s}^2}}=0$ (Eyov et al. Reference Eyov, Klar, Kadri and Stiassnie2013). The first progressive mode ($n=1$) for an elastic seabed is a Scholte wave, which propagates all the way to the shore, where it turns into a Rayleigh wave. From (30) of Eyov et al. (Reference Eyov, Klar, Kadri and Stiassnie2013), the critical frequency for a particular depth is given by

(4.1)\begin{equation} \omega_{en}=\left(n-\frac{3}{2}\right){\rm \pi} \frac{C_{l}C_{s}}{h_{en} \sqrt{C_{s}^2 - C_{l}^2}}\quad n=2,3,\ldots, \end{equation}

which is a good approximation to the actual critical frequency, though it is based on the location of the vertical asymptotes in the dispersion relation plot – location of red circle in figure 6(a). For accuracy, we require a better approximation to the actual intersection of the two curves in the dispersion relation plot – blue circle in figure 6(a). We begin with the dimensionless form of the dispersion relation (3.44):

(4.2)\begin{equation} \tanh(\tilde{r})=\frac{\dfrac{\tilde{\omega}^{2}}{\tilde{r}} (\varepsilon_{2}+\varepsilon_{1})}{\dfrac{\tilde{\omega}^{4}}{\tilde{r}^{2}} \varepsilon_{2}+\varepsilon_{1}},\quad \varepsilon_{1}=\frac{4\tilde{k}^{2}\tilde{q}\tilde{s} \tilde{\mu}}{\tilde{k}^{2}+\tilde{s}^{2}}+\tilde{\lambda} \tilde{k}^{2}-(\tilde{\lambda}+2\tilde{\mu})\tilde{q}^{2}, \quad \varepsilon_{2}=\tilde{q}\left(\frac{\tilde{k}^{2}- \tilde{s}^{2}}{\tilde{k}^{2}+\tilde{s}^{2}}\right), \end{equation}

with $r, k, q, s, \omega, \lambda, \mu, C_{l}, C_{s}, C_{p}$ made dimensionless according to

(4.3)\begin{equation} \left.\begin{array}{c} \displaystyle \tilde{r}=hr, \quad \tilde{k}=hk, \quad \tilde{q}=hq, \quad \tilde{s}=hs, \quad \tilde{\omega}=\sqrt{\dfrac{h}{g}}\omega , \quad \tilde{\lambda}=\dfrac{\lambda}{\rho_{l} gh}, \quad \tilde{\mu}=\dfrac{\mu}{\rho_{l} gh}, \\ \displaystyle \widetilde{C_{l}}=\dfrac{C_{l}}{\sqrt{gh}}, \quad \widetilde{C_{s}}=\dfrac{C_{s}}{\sqrt{gh}}, \quad \widetilde{C_{p}}=\dfrac{C_{p}}{\sqrt{gh}}. \end{array}\right\} \end{equation}

We substitute

(4.4a–c) \begin{equation} \tilde{k}=\sqrt{\tilde{r}^{2}+\frac{\tilde{\omega}^{2}}{\widetilde{C_{l}}^{2}}}, \quad \tilde{q}=\sqrt{\tilde{r}^{2}+\frac{\tilde{\omega}^{2}}{\widetilde{C_{l}}^{2}}- \smash{\frac{\tilde{\omega}^{2}}{\widetilde{C_{p}}^{2}}}}, \quad \tilde{s}=\sqrt{\tilde{r}^{2}+\frac{\tilde{\omega}^{2}}{\widetilde{C_{l}}^{2}}- \frac{\tilde{\omega}^{2}}{\widetilde{C_{s}}^{2}}} \end{equation}

into (4.2) to obtain a function of $\tilde {r}$ alone, followed by another substitution $\tilde {r}\Rightarrow \textrm {i} \tilde {r}$ to retrieve the acoustic–gravity wave solutions.

Let $\tilde {\omega }={\tilde {r} \widetilde {C_{s}} \widetilde {C_{l}}}/{\sqrt {\widetilde {C_{s}}^{2}-\widetilde {C_{l}}^{2}}}$, which is the $\tilde {s}=0$ condition for the termination of progressive modes, and then let $\tilde {r}=(n-\frac {3}{2}){\rm \pi} +\delta (n)$, where $\delta (n)$ represents a small (mode-dependent) positive offset away from the vertical asymptotes located at $\tilde {r}=(n-\frac {3}{2}){\rm \pi}$. So the desired $\delta (n)$ is the $\tilde {r}$ separation between the blue and red circles in figure 6(a). Ignoring terms of $O(\delta (n)^2)$ an approximation of the dispersion relation (see Appendix B) can be written as

(4.5)\begin{equation} -\cot(\delta(n))=\frac{a_{n}+b_{n}\delta(n)}{c_{n}+d_{n}\delta(n)}, \end{equation}

where the coefficients $a_{n}, b_{n}, c_{n}, d_{n}$ are given in Appendix C. Then using the approximation $-\cot (\delta (n))\simeq -{1}/{\delta (n)}$ for small $\delta (n)$, (4.5) can be put into quadratic form:

(4.6)\begin{equation} b_{n}\delta(n)^{2}+(a_{n}+d_{n})\delta(n)+c_{n}=0. \end{equation}

This can be solved for $\delta (n)$, and then the value of $\tilde {r}_{n}$ and the critical frequency $\omega _{n}$ can be obtained from

(4.7)\begin{equation} \tilde{r}_{n}=\left(n-\frac{3}{2}\right){\rm \pi}+\delta(n). \end{equation}

To determine how well $\delta (n)$ predicts the offset a comparison was made between the approximate value of $\omega _{n}$ calculated using (4.7) and that found by using the dispersion relation

(4.8)\begin{equation} {\rm error [\%]}=\left|\frac{\omega_{n}({\rm approx})- \omega_{n}({\rm dispersion})}{\omega_{n}({\rm dispersion})} \times100\right|. \end{equation}

The comparison was carried out for depths ranging from 500 to 8000 m and all available modes. The maximum error occurred in the second mode at a depth of 8000 m, but was still less than 0.1 % (see figure 7). In § 5 the results for $\tilde {r}_{n}$ and $\delta (n)$ are used to construct approximate phase velocity curves.

Figure 7. Percentage error for approximate critical frequencies $\omega _{n}$ from (4.7). Depths range between 500 m (lower error bound) and 8000 m (upper error bound) – all available modes.

4.2. Surface wave

The surface gravity mode 00 does not exist for all frequencies in the elastic seabed case and never exists for a rigid seabed. The cut-off condition for mode 00 is given by (3.66) which can be solved numerically for a given depth. Alternatively a good approximation can be obtained by seeking the frequency at which the gradient of the left-hand side of the dimensionless dispersion relation in (4.2) is equal to the gradient of the right-hand side. This occurs for small (ultimately zero) $\tilde {r}$. In this case we make the approximation $\tanh (\tilde {r})\simeq \tilde {r}$. Now we differentiate (4.2) with respect to $\tilde {r}$. Then we express the result as a series in $\tilde {r}^2$ to arrive at

(4.9)\begin{equation} 1=\frac{1}{\tilde{\omega}^2}+\frac{\mathcal{A}}{\tilde{\omega}} +O(\tilde{r}^2). \end{equation}

In the limit $\tilde {r}\rightarrow 0$, (4.9) can be written as the quadratic $\tilde {\omega }^2 -\mathcal {A} \tilde {\omega } -1=0$, with

(4.10)\begin{align} \mathcal{A}=2\frac{\left[\left(\tilde{\mu}+\dfrac{\tilde{\lambda}}{2}\right) \tilde{C}_{l}^2-\tilde{\mu}\tilde{C}_{p}^2\right] (2\tilde{C}_{s}^2-\tilde{C}_{l}^2)\sqrt{\tilde{C}_{s}^2- \tilde{C}_{l}^2}+2\tilde{C}_{p}\tilde{C}_{s}\tilde{\mu} \sqrt{\tilde{C}_{p}^2-\tilde{C}_{l}^2} (\tilde{C}_{s}^2-\tilde{C}_{l}^2)}{\tilde{C}_{l}^3 \tilde{C}_{p} \sqrt{\tilde{C}_{p}^2-\tilde{C}_{l}^2} \sqrt{\tilde{C}_{s}^2-\tilde{C}_{l}^2}}. \end{align}

Taking the positive root of the quadratic gives the approximate cut-off frequency which we name $\tilde {\varOmega }_{00}$. The dimensional form can be recovered from $\varOmega _{00}=\tilde {\varOmega }_{00}\sqrt {\vphantom{\frac{a}{a}}{g}/{h}}$. A workable approximation to the cut-off frequency can be obtained by taking $\mathcal {A}$. We call this approximation $\mathcal {A}_{00}$. Table 1 gives values for $\omega _{00}$ found using a numeric solver and compares with the approximations $\varOmega _{00}$ and $\mathcal {A}_{00}$ indicating errors for various depths. The error values were calculated from

(4.11)\begin{equation} {\rm error [\%]}=\left|\frac{\varOmega_{00}({\rm or } \mathcal{A}_{00})-\omega_{00}}{\omega_{00}}\times100\right|. \end{equation}

Table 1. Comparison of cut-off frequencies obtained from numeric solver ($\omega _{00}$) with approximations from quadratic solution ($\varOmega _{00}$) and coarse approximation ($\mathcal {A}_{00}$) for various depths $h$.

6.1. Acoustic–gravity waves

Figure 13 compares the first acoustic mode for the elastic case (3.92) and rigid case ((3.22) of Williams et al. (Reference Williams, Kadri and Abdolali2021)). As in Eyov et al. (Reference Eyov, Klar, Kadri and Stiassnie2013) the values used for $\rho _{l}, \rho _{s}, C_{l}, C_{s}$ and $C_{p}$ are average values taken from Dziewonshi & Anderson (Reference Dziewonshi and Anderson1981). One stark difference between the two cases is that the signal terminates after some time in the elastic case, whereas it continues indefinitely in the rigid case. Another difference is the presence of signal at times earlier than the main pulse in the elastic case, but no signal at all in the rigid stationary phase model. The phase velocity curves for the elastic case (figure 11) indicate that frequencies close to the critical frequency for each mode receive a boost in phase velocity enabling signals to propagate faster. For these frequencies speeds close to $C_{s}$ are achievable. The rigid seabed stationary phase model produces complex numbers for times earlier than ${x}/{C_{l}}$ due to a singularity induced by the stationary phase method (Stiassnie Reference Stiassnie2010). The pressure amplitudes are similar.

Figure 13. First acoustic mode (a) with elastic seabed and (b) with rigid seabed.

A sensitivity analysis was carried out looking at the effects of six parameters on the signal duration (see figure 14). Each parameter was varied individually away from its reference value (table 2) while holding all other parameters at reference. Then the percentage change in pulse duration was divided by the percentage change in the parameter to arrive at the sensitivity value. It was found that the rigidity of the seabed most affected the signal duration. Increasing the Lamé parameters increases $C_{s}$ and $C_{p}$ in accordance with (2.7a,b); the ratio ${C_{p}}/{C_{s}}$ was kept constant. Another difference between the elastic seabed case and the rigid seabed case appears when a fast Fourier transform of the signals is examined. The elastic case shows a slight upward shift in the frequency peak (see figure 15). This is in contrast to the slight downward shift found when a viscous compressible sediment layer is overlying the seabed as in Abdolali, Kirby & Bellotti (Reference Abdolali, Kirby and Bellotti2015). Whilst investigating the synthetic acoustic–gravity waves, we found that band-pass filtering applied to the signal generated by combining the first 10 modes (figure 16a) revealed some interesting peaks located close to the expected arrival time for a phase velocity of $C_{l}$ (figure 16b). There are four peaks of particular interest labelled 1, 2, 3 and 4. The presence of the peaks is a consequence of the fault's geometry and motion, and is not related to the rigidity of the seabed, since the peaks are also present under rigid seabed conditions, and even when the signal considered was purely acoustic, as in Mei & Kadri (Reference Mei and Kadri2018). The time spacing between pairs of peaks responds to changes in either fault half-width $b$ or rupture duration $\tau =2T$ in a linear fashion, so that details of the fault's geometry and dynamics are encoded in the acoustic–gravity waves. Time $\Delta t_{1}$ between peak numbers 1 and 2 (or 3 and 4) is exactly the rupture duration, and $\Delta t_{2}$ between peaks 1 and 3 (or 2 and 4) is linearly related to the fault half-width through $\Delta t_{2}={2b}/{C_{l}}$. When the slender fault begins to move, peaks 1 and 3 are generated at the edges of the slender fault and begin to propagate. The time separation between these peaks is explained as the time required for a wave travelling at speed $C_{l}$ to cross a fault width of $2b$. At the end of the fault's motion, after $\tau$ seconds, the second pair of peaks (2 and 4) is generated and propagates away – also separated in time by $\Delta t_{2}$. The resultant waveform as would be seen in the far field is a collection of four peaks. The amplitude of the peaks depends linearly on the uplift velocity $W_{0}$. The timings between peaks agree well with the values given in table 2. Let subscripts $1,2,3$ and $4$ represent the peaks denoted by the numbers 1, 2, 3 and 4 respectively. Then, $\Delta t_{12}=10.97\ \textrm {s}$, $\Delta t_{34}=8.9\ \textrm {s}$, $\Delta t_{13}=57.67\ \textrm {s}$ and $\Delta t_{24}=55.46\ \textrm {s}$. The times $\Delta t_{12}$ and $\Delta t_{34}$ represent $\tau$, the uplift time, which (from table 2) is actually 10 s. The times $\Delta t_{13}$ and $\Delta t_{24}$ are the transit times for an acoustic signal to cross a fault width $2b$, which, again from table 2, is actually 55.17 s. The information that could be extracted from the timings embedded in the acoustic modes would be of interest to the inverse process that reconstructs fault parameters from received signals (Gomez & Kadri Reference Gomez and Kadri2021). In § 6.3 an actual hydrophone recording made during the Samoa 2009 event is filtered to reveal the four peaks encoded within it. Nosov (Reference Nosov.1999) concludes that the acoustic modes have a frequency spectrum which depends on the time history and spatial structure of the bottom displacement – which is referred to as the tsunami's voice. This is exactly what we find encoded into the characteristic (four) peaks.

Figure 14. Response of signal duration when changing parameters.

Figure 15. Fast Fourier transform of first four available modes. Depth $h= 4000$ m.

Figure 16. Band-pass filtering applied to the 10 combined modes of the synthetic acoustic–gravity wave generated by a single slender fault. (a) The first 10 modes combined. (b) The resulting signal after application of band-pass filtering. The characteristic peaks are numbered 1, 2, 3 and 4.

Table 2. Constants and parameters used in comparison of elastic seabed with rigid seabed.

6.2. Surface waves

Consider now the surface waves generated by the single slender fault with parameters as per table 2. The equations generating the surface waves are (3.93) for the elastic case and (3.23) of Williams et al. (Reference Williams, Kadri and Abdolali2021) for the rigid case. At a depth of $h=4000$ m there is little difference between elastic and rigid cases (figure 17a), but at a depth of $h=1000$ m differences are more apparent (figure 17b). In deeper water the surface wave is almost unaffected by the elasticity of the seabed (Eyov Reference Eyov2013). This surface wave is the main tsunami, i.e. mode 01.

Figure 17. Surface elevation comparison (elastic versus rigid). Coordinates are $x=1000$ km, $y=0$ km. Coordinate origin at fault centroid. Depths (a) $h=4000$ m and (b) $h=1000$ m.

When the seabed is elastic the possibility of a second surface wave arises. This wave does not exist for all frequencies and never exists in the rigid case (Eyov Reference Eyov2013). The gradient condition (3.66) has to be satisfied before mode 00 can propagate (see figure 18). The mode 00 surface wave has a phase velocity $C_{l}$, and a negligible amplitude, of the order of micrometres. A plot of mode 00 under the conditions of table 2 can be found in figure 19. In the plot there are four distinct peaks numbered 1 to 4. These peaks in the mode 00 surface wave correspond to the peaks numbered 1, 2, 3 and 4 found in the acoustic signal. The assumption of a rectangular fault moving at a uniform speed results in symmetric peaks, whereas in reality the motion is much more complicated thus upsetting the symmetry of the peaks seen in figure 19.

Figure 18. Left-hand side (LHS) of dispersion relation (3.44) (dashed trace) and right-hand side (RHS) of dispersion relation (solid trace) when $r\in \mathbb {R}$. The frequency is at the point where mode 00 becomes active, $\omega \simeq 6.95\ \textrm {rad s}^{-1}$, $h=4000$ m (see table 1). The top logarithmic plot indicates overall behaviour. The middle and bottom plots provide expanded views. The mode 00 solution in the middle plot is where the solid curve touches the dashed curve $0\leq r\leq 0.1$, and the mode 01 solution (the usual tsunami) in the bottom plot is where the solid curve again makes contact with the dashed curve in the descending phase $2.5 \leq r \leq 3$.

Figure 19. Mode 00 surface–gravity wave with envelope.

6.3. Hydrophone recordings

The theory developed in this paper leading to the equations for pressure (3.92) and surface elevation (3.93) is linear. Therefore, as in Williams et al. (Reference Williams, Kadri and Abdolali2021), more complicated multi-fault scenarios can be constructed from single slender fault solutions by linear superposition, given that the parameters for each individual fault are known. To contrast the case of an elastic seabed with that of a rigid seabed we revisit the Sumatra 2004 earthquake discussed in § 5.1.2 of Williams et al. (Reference Williams, Kadri and Abdolali2021). The geographical area considered here ranges over [$70^{\circ }$ to $100^{\circ }$] E longitude and [$-15^{\circ }$ to $20^{\circ }$] N latitude. This curved patch on the (idealised) spherical Earth is mapped to flat $x, y$ coordinates. The conversion factor of metres per degree is fixed for the latitude ($y$) direction, but the metres per degree in the longitudinal ($x$) direction varies with latitude, being maximum at the equator and decreasing as the poles are approached. To simplify the calculations an average value for metres per degree longitude was used and the area considered was kept reasonably small.

The fault centroids are marked by black stars in figures 20(a) and 20(b) and all faults are contained within the masked off ‘earthquake zone’. The purpose of the earthquake zone was to avoid pressure calculations too close to the faults. The location of hydrophone H08N is marked with a red star. Figure 20 indicates the time evolution of the bottom pressure signal for both rigid and elastic seabeds. The elastic seabed has pressure signals already close to the hydrophone at $t=1000$ s, whereas the rigid seabed only has pressure signals local to the earthquake zone at this time. This is due to the twin effects of a boost in phase velocity for frequencies close to critical in the elastic case and the absence of signal ahead of the main pulse in the rigid stationary phase model. As time proceeds the pressure signals for the elastic case can be seen to traverse the area considered, so that by $t=3625$ s the area is largely clear of pressure oscillations. In contrast, the rigid stationary phase model shows persistent and ever-increasing pressure oscillations around the earthquake zone. The elastic seabed could therefore be considered more physically realistic.

Figure 20. Bottom pressure comparison between rigid and elastic seabed. The location of H08N hydrophone is indicated by a red star at bottom left of each panel. By 3625 s, the elastic model has largely cleared of acoustic–gravity waves whereas the rigid model still has strong oscillations around the earthquake zone. (a) Rigid seabed, bottom pressure map calculated at nine time intervals after first fault movement. (b) Elastic seabed, bottom pressure map calculated at the same time intervals.

Figure 21 compares the predictions made by the elastic seabed model against recorded data for the Sumatra event derived from the southern (H08S1) hydrophone and the seismograph at nearby Diego Garcia (DGAR). The signals recorded by the three hydrophones at station H08S were very similar to each other, and so only that of H08S1 is displayed in the plots (similarly for station H08N). The amplitude of the main acoustic–gravity wave signal for the northern triad is much smaller than that of the southern triad – possibly due to the shielding effect of the Chagos Archipelago (see figures 22a and 23). However, the leading pulses ($P$-waves) are of similar amplitude. This suggests that the detection of the $P$-waves by the hydrophones is largely unaffected by the presence of the island – unlike signals that travel only through the water. The hydrophone data were obtained from the Comprehensive Nuclear Test Ban Treaty Organisation (CTBTO) and the seismic data from the Incorporated Research Institutions for Seismology (IRIS). The fault configuration is that of table 3/figure 6(a) (Williams et al. Reference Williams, Kadri and Abdolali2021) and the start time given by the United States Geological Survey (USGS) website is UTC 2004-12-26 00:58:53, which corresponds to $t=0$. In the plots the blue vertical lines correspond to the expected arrival time of acoustic–gravity waves travelling at a phase speed of $1450\ \textrm {m s}^{-1}$. The green vertical lines correspond to acoustic–gravity waves travelling at a phase speed of $C_{s}=3550\ \textrm {m s}^{-1}$ and the red vertical lines correspond to a phase speed of $\simeq 8000\ \textrm {m s}^{-1}$. The DGAR seismometer records small-amplitude $P$-waves arriving at the red vertical line, which then transition to larger-amplitude $S$-waves at the green line. The hydrophone H08S weakly detects the $P$-wave activity. The seismic $S$-waves do not exist in the liquid (no shear). The main acoustic–gravity wave signal then arrives and is detected by the hydrophone at the blue line. The re-scaled plot of the hydrophone signal shows this behaviour more clearly. Since $C_{s}$ is the speed limit for the acoustic–gravity waves in the elastic model, the model is unable to predict the $P$-wave portion of the hydrophone signal. Modifications to the existing model or maybe a new model would be needed to capture this behaviour. Between the green and blue lines the elastic model predicts acoustic–gravity waves that can travel at phase speeds close to $C_{s}$ for some frequencies. The hydrophones show weak signal in this region, possibly due to the filtering effect of the hydrophone's response. After the blue line the elastic model predicts a signal that is close in amplitude to that of the hydrophone recordings, but decays more slowly – possibly due to a lack of dissipation included in the model. However, the signal duration is at least finite in the elastic case. There are processes missing from the elastic model (varying bathymetry, reflection, refraction, dissipation etc.) so an exact match is not expected. Examination of figure 21 shows the arrival times for $P$-waves, $S$-waves and the main acoustic–gravity wave pulse (travelling at a phase speed of $C_{l}$) are consistent with our assumptions of constant water density, constant speed of sound in liquid and constant speed of propagation in solid. The leading pulse seen in the hydrophone recordings is primarily made up of lower-frequency components. To show this, a band-pass filter was applied to the hydrophone recording at H08S. The filter eliminates most of the frequency components below 3 Hz and has largely flattened the leading pulse of the H08S signal (figure 24). In order to enhance the detection of acoustic signal between the red and blue lines, ultralow-frequency hydrophones should be used.

Figure 21. Comparison of the current elastic model with both hydrophone and seismic data for the Sumatra 2004 event. The time axis begins at UTC 2004-12-26 00:58:53 ($t=0$). The vertical red line represents the arrival time for a propagation speed of $8000\ \textrm {m s}^{-1}$, the vertical green line represents the arrival time for a propagation speed of $C_{s} = 3550\ \textrm {m s}^{-1}$ and the vertical blue line represents the arrival time for a propagation speed of $C_{l}=1450 \ \textrm {m s}^{-1}$.

Figure 22. (a) Locations for the H08N and H08S hydrophone triads, along with the DGAR seismograph (yellow markers). The northern triad is shielded by the Chagos Archipelago. (b) Expanded view of island and west coast of Sumatra. Images from Google Earth.

Figure 23. (a,b) Overlay of elastic model prediction onto hydrophone data of north and south locations. (c,d) North and south hydrophone data with re-scaled vertical axis. Red vertical line, arrival time for phase speed $8000\ \textrm {m s}^{-1}$; green vertical line, arrival time for phase speed $C_{s} = 3550\ \textrm {m s}^{-1}$; blue vertical line, arrival time for phase speed $C_{l} = 1450\ \textrm {m s}^{-1}$.

Figure 24. Leading pulse of hydrophone signal is largely made up of low-frequency components which filtering is able to suppress.

Table 3. Comparison of two key fault parameters (rupture duration and width) obtained by different methods. The second column ($\Delta t_{1,2,3,4}$) reports data obtained by filtering the H11 hydrophone signal and measuring timings between peaks. The third column reports data obtained by the methods described within Gomez (Reference Gomez2022). The data in the fourth column are estimates derived from the USGS website.

To demonstrate the extraction of fault timing and geometry from acoustic–gravity wave signals a band-pass filter was applied to the data obtained from hydrophone H11 located at Wake Island during the Samoa 2009 event (data supplied by CTBTO). The timings for the peaks are $\Delta t_{12}=15.46$, $\Delta t_{34}=27.75$, $\Delta t_{13}=29.89$ and $\Delta t_{24}=42.18\,\textrm {s}$ (figure 25). Note that the time axis in figure 25 does not begin from the start of the rupture. The time axis in this case represents an 1800 s window around the main hydrophone signal. This is not a concern here since only time differences ($\Delta t$ values) are required. Unlike in the synthetic case this event is asymmetric in the sense that the timings suggest either a trapezoidal rupture geometry or non-uniform uplift velocity (or both). The time $\Delta t_{12}=15.46$ s represents the uplift time for the leading edge of the fault. Assuming that the front and back edges of the fault begin moving together, then $\Delta t_{13}=29.89$ s represents the time for the acoustic signal to travel from the back of the fault, across the fault width, to the front and thus indicates a fault width of $2b=C_{l}\Delta t_{13}\equiv 43.3$ km. The time $\Delta t_{34}=27.75$ s represents the total time for the fault movement (the back end continues moving after the front has stopped). These data compare quite well with those retrieved via inverse modelling in Gomez (Reference Gomez2022). Also, the fault width and timing are approximately those found in the USGS finite fault model (see figure 26) at https://earthquake.usgs.gov/earthquakes/eventpage/usp000h1ys/executive.

Figure 25. (a) Recorded hydrophone data from H11 at Wake Island for Samoa 2009 event. Note that $t=0$ does not correspond to the rupture start time. (b) Signal after application of band-pass filtering, focusing on the time interval containing the initiation of the main pulse. Data sampling occurs at 250 Hz (one sample every 4 ms).

Figure 26. The USGS finite fault model dimensions and timings. (a) The USGS finite fault model for Samoa 2009 event with scale at bottom left corner. Rectangular region shown is approximately $30\ \textrm {km} \times 180\ \textrm {km}$. (b) The USGS moment rate function for Samoa 2009 event. The main peak ends between 25 and 35 s.

6.4. The DART buoy data

For the validation of the surface wave calculations against real data we consider the Tohoku event of March 2011 as covered in Williams et al. (Reference Williams, Kadri and Abdolali2021). The parameters used in the elastic model of this paper were changed slightly from those found in Williams et al. (Reference Williams, Kadri and Abdolali2021) and are listed in table 4. In Williams et al. (Reference Williams, Kadri and Abdolali2021) the event was treated as a multi-fault event, so as to capture the main tsunami. However, this paper uses the elastic model and the middle term of (3.93), integrated directly to describe the tsunami (mode 01). It was not necessary to split the fault into a number of faults. By integrating directly the tsunami could be modelled by a single fault. The main peak of the tsunami is described quite well by the elastic model, in terms of both timing and amplitude (see figure 27).

Table 4. Constants and parameters used in the calculation of surface elevation at DART buoy 21418 for Tohoku 2011 event – elastic model. Also refer to Williams et al. (Reference Williams, Kadri and Abdolali2021).

Figure 27. Surface elevations compared for Tohoku 2011 event at DART buoy 21418.