2. Waves in the ocean–seafloor system

Waves in the two media are coupled together by the boundary conditions at the interface. The two media are considered to be mildly non-uniform, in that the acoustic speed in ocean water is a slowly varying function of the vertical coordinate, and the acoustic speeds in the seafloor are slowly varying functions of the horizontal coordinates (see supplementary material § 1.1). For the present analysis, the water depth $h$ is assumed to be uniform. Although a typical seafloor may be comprised of several layers that could play a role in the overall dynamics, here the seafloor is assumed to be made up of a single semi-infinite layer with averaged properties. Energy dissipation effects in water and along the seafloor are assumed to be negligible. In this and the following sections, consistent with the representations in Phillips (Reference Phillips1960) and Hasselmann (Reference Hasselmann1963), the symbols ${\rm d}C_1$, ${\rm d}C_2$, ${\rm d}A_1$, ${\rm d}A_2$, ${\rm d}A'_1$, ${\rm d}B'_1$ are used to denote the Fourier coefficients in the Fourier–Stieltjes integral representations of the stationary random waves being analysed.

As shown in figure 1, the origin is located on the seafloor, with $(x_1, x_2)$ representing two mutually orthogonal horizontal directions, and with the $y$ axis pointing vertically upwards. However, a more insightful and compact representation is possible if a particular wave group is followed, with the $x$ axis along the propagation direction as represented by a wavenumber ${\boldsymbol k}$. For a propagation direction $\theta _p$ in the $(x_1, x_2)$ plane,

(2.1a) \begin{gather} x = x_1\cos\theta_p + x_2\sin\theta_p, \end{gather}

(2.1b) \begin{gather} k_1 = k\cos\theta_p,\quad k_2 = k\sin\theta_p, \end{gather}

(2.1c) \begin{gather} {\rm and}\quad {\boldsymbol k}\boldsymbol{\cdot} {\boldsymbol x} = kx,\quad k^2 = k^2_1 + k^2_2. \end{gather}

Equation (2.1b) follows because the wave group has no propagation in the direction perpendicular to $\theta _p$.

Figure 1.

Schematic sketch (not to scale) providing an overview of the system being studied, along with some definitions.

The acoustic-gravity–Scholte wave system is constrained by boundary conditions at $y = 0$ and $y = h$, and the radiation conditions at $x\rightarrow \pm \infty$. Here, the attention is on an acoustic-gravity–Scholte wave system propagating from left to right along the positive $x$ direction. Hence for $x \rightarrow -\infty$, the waves are incoming, while as $x\rightarrow \infty$, the radiation condition for outgoing waves applies. Viscous forces in water are neglected in this analysis (Kibblewhite Reference Kibblewhite1989), and the water column is thought to support only compression/extension waves but not shear waves (Webb Reference Webb1986; Kibblewhite Reference Kibblewhite1989). Energy dissipation over the seafloor is also neglected (but see Kibblewhite Reference Kibblewhite1989).

The vertical component of the particle velocity at $y = 0$ must be equal across the boundary, so that water particles adjacent to the seafloor move with the same vertical velocity as the solid particles attached to the seafloor. In addition, the normal stresses on the seafloor on the water side and the seafloor side must be equal. The seafloor supports both extensional (dilatational) waves and shear (distortional) waves. Hence since the water column does not support shear waves, the shear stress in the solid at the seafloor must be zero. At the top of the water column, the pressure under quiescent conditions just equals the atmospheric pressure. Since the seafloor also feels this atmospheric pressure, it can be subtracted out of the normal stress on both $y= 0$ and $y=h$ surfaces. Compression waves through the water column and seafloor waves can exist even in the absence of excitation from below or above. These waves are referred to as free waves, and may owe their existence to impulsive forcing due to explosions or past surface-wave activity. Their properties are constrained by the properties of the two media, their frequencies and wavenumbers being related together by a dispersion relation. The immediate objective below is to derive the dispersion relation for the two-media system that is consistent with the interface boundary conditions.

The treatment below is based on Graff (Reference Graff1991, chapter 6). Compression wave motion in the water column is assumed to be irrotational and hence is described by a scalar displacement potential $\phi (x,y;t)$. On the other hand, the motion in the solid is not purely irrotational and in general includes a rotational component because the solid can support shear deformation. Hence the solid medium requires a scalar potential $\phi '(x,y;t)$ and an additional vector potential ${\boldsymbol H}'(x,y,z;t)$. The purpose here is to study a wave system propagating in the $(x,y)$ plane, so ${\boldsymbol H}'$ is restricted to be ${\boldsymbol H}' = (0, 0, H'_z)$.

The displacement ${\boldsymbol u}(x,y;t)$ as represented by components $(u,v)$ of the water particles is given by

(2.2a,b)\begin{equation} {\boldsymbol u} = \boldsymbol{\nabla} \phi\quad \Rightarrow\quad u = \frac{\partial \phi}{\partial x}\ {\rm and}\ v = \frac{\partial \phi}{\partial y}. \end{equation}

The displacement ${\boldsymbol u}'(x,y;t)$ in the solid is given by

(2.3)\begin{equation} {\boldsymbol u}' = \boldsymbol{\nabla} \phi' + {\boldsymbol \nabla}\times {\boldsymbol H}'. \end{equation}

The horizontal and vertical components of ${\boldsymbol u}'$ are given by

(2.4a,b)\begin{equation} u' = \frac{\partial \phi'}{\partial x} + \frac{\partial H_z}{\partial y}\quad {\rm and}\quad v' = \frac{\partial \phi'}{\partial y} - \frac{\partial H_z}{\partial x}. \end{equation}

The kinematic boundary condition at the seafloor is

(2.5)\begin{equation} v = v'\quad \Rightarrow\quad \frac{\partial \phi}{\partial y} = \frac{\partial \phi'}{\partial y} - \frac{\partial H_z}{\partial x},\quad {\rm at}\ y = 0. \end{equation}

The stress component due to the acoustic-gravity waves in the water that is normal to the seafloor is

(2.6)\begin{equation} \tau_{yy} = \lambda\left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right), \end{equation}

while the equal and opposing normal stress component in the solid is (cf. Eyov et al. Reference Eyov, Klar, Kadri and Stiassnie2013)

(2.7)\begin{equation} \tau'_{yy} = \lambda'\left(\frac{\partial u'}{\partial x} + \frac{\partial v'}{\partial y}\right) + 2\mu'\,\frac{\partial v'}{\partial y}. \end{equation}

Here, $\lambda$ denotes the Lamé constant for water, and $\lambda '$ and $\mu '$ represent the Lamé constants for the seafloor; $\lambda$, $\lambda '$ and $\mu '$ may all be slowly varying spatially, consistent with the discussion in supplementary material § 1.1. The first of the kinetic boundary conditions at the seafloor is

(2.8)\begin{equation} \tau_{yy} = \tau'_{yy},\quad {\rm at}\ y = 0. \end{equation}

The shear stress component along the propagation direction in the $y$-plane on the seafloor is given by

(2.9)\begin{equation} \tau'_{xy} = \mu'\left(\frac{\partial u'}{\partial y} + \frac{\partial v'}{\partial x}\right). \end{equation}

The second kinetic boundary condition on the seafloor is

(2.10)\begin{equation} \tau'_{yx} = \tau'_{xy} = \mu'\left(\frac{\partial u'}{\partial y} + \frac{\partial v'}{\partial x}\right) = 0,\quad {\rm at}\ y = 0. \end{equation}

A kinetic boundary condition can also be specified for the top of the water column, $y = h$. For free-wave propagation, the pressure at $y= h$ is set to zero, with the static atmospheric pressure subtracted out. For free waves, in terms of the components $(u,v)$, this condition can be expressed as

(2.11)\begin{equation} \lambda \left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right) = 0,\quad y = h. \end{equation}

As mentioned, the free-surface-wave field causing the present acoustic-gravity–Scholte waves can be described as a stationary random process (Hasselmann Reference Hasselmann1963). With the acoustic-gravity–Scholte wave system also treated as a stationary random process, a Fourier–Stieltjes integral representation is used here to expand the potentials $\phi$, $\phi '$ and ${\boldsymbol H}'$ (Phillips Reference Phillips1977, chapter 4). The acoustic-gravity wave field in the water column for rightward propagating waves can be expressed as

(2.12)\begin{equation} \phi(x,y;t) = \iint_{k,\omega}\bigl({\rm d}{C_1 \exp(}{{\rm i}}\alpha y) + {\rm d}C_2\exp({-{\rm i}\alpha y})\bigr)\exp({{\rm i} (k x - \omega t)}). \end{equation}

It is recalled that $kx = {\boldsymbol k}\boldsymbol {\cdot }{\boldsymbol x}$. The acoustic-gravity waves travel along inclined paths, with $\textrm {d}C_1 = \textrm {d}C_1(\omega,{\boldsymbol k})$ and $\textrm {d}C_2 = \textrm {d}C_2(\omega, {\boldsymbol k})$ representing, respectively, the Fourier components of the potential defining the acoustic-gravity waves propagating towards the seafloor, and the Fourier components for the acoustic-gravity waves reflected from the seafloor. These waves propagate as pressure waves, and in a broad sense can be considered as guided waves. Here, $\omega$ is the frequency of the acoustic-gravity and Scholte waves. When these waves result from a second-order interaction of surface waves of frequencies $\sigma _1$ and $\sigma _2$, we have $\omega = \sigma _1 + \sigma _2$. Here, ${\boldsymbol k}$ denotes the horizontal wavenumber, and $\alpha$ is the vertical wavenumber, given by

(2.13)\begin{equation} \alpha = \sqrt{(\omega^2/c^2_1 - k^2)}, \end{equation}

where $c_1$ represents the phase velocity of acoustic waves in the water column, and is taken to be a mildly varying function of the depth coordinate $y$. Supplementary material § 1.1 discusses the type of variations for which such an approximation would be admissible, and considers examples where the vertical distribution of acoustic speed in the ocean is described by the Munk profile (Munk Reference Munk1974).

The combined wave field (acoustic-gravity and Scholte) is free to propagate in the horizontal direction, and here the radiation conditions consistent with left-to-right propagation are used. In the water column, the acoustic-gravity wave field is trapped vertically between two boundaries, $y= 0$ and $y = h$. Hence $\phi (x,y;t)$ can be represented alternatively as

(2.14)\begin{equation} \phi(x,y;t) = \iint_{k,\omega} ({\rm d}A_1 \cos\alpha y + {\rm d}A_2 \sin\alpha y)\exp({{\rm i} (k x - \omega t)}), \end{equation}

where $\textrm {d}A_1 = (\textrm {d}C_1 + \textrm {d}C_2)$ and $\textrm {d}A_2 = \textrm {i}(\textrm {d}C_1 - \textrm {d}C_2)$. Both representations are used in this work.

In general, the wave field in the seafloor could be comprised of bulk waves and surface waves. Bulk-wave components possible in the present situation are one seismic pressure wave (P-wave) component with particle oscillation along the propagation direction, and one shear wave (SV-wave) component with particle oscillations perpendicular to the propagation direction but in the $(x,y)$ plane (Graff Reference Graff1991). Since the seafloor layer extends to $y \rightarrow -\infty$, the Fourier–Stieltjes representations for the potentials $\phi '$ and $H'_z$ in the solid are

(2.15)\begin{equation} \phi'(x,y; t) = \iint_{k, \omega} {\rm d}A'_1 \exp({-{\rm i}\bar{\alpha} y}) \exp({{\rm i}(kx - \omega t)}) \end{equation}

and

(2.16)\begin{equation} H'_z(x,y; t) = \iint_{k, \omega} {\rm d}B'_1 \exp({-{\rm i}\bar{\beta}y})\exp({{\rm i}(kx - \omega t)}), \end{equation}

where $\textrm {d}A'_1 = \textrm {d}A'_1(\omega,{\boldsymbol k})$ and $\textrm {d}B'_1 = \textrm {d}B'_1(\omega,{\boldsymbol k})$ are the Fourier coefficients representing the extensional (dilatational) or P-wave and distortional SV-wave potential components, respectively, through the seafloor (Graff Reference Graff1991, chapter 6). Thus

(2.17a,b)\begin{equation} {\bar{\alpha}} = \sqrt{(\omega^2/c^2_p - k^2)}\quad {\rm and}\quad {\bar{\beta}} = \sqrt{(\omega^2/c^2_s - k^2)}, \end{equation}

where $c_p$ and $c_s$ denote, respectively, the phase velocities of extensional and shear waves in the seafloor, and are taken to be mildly varying functions of the horizontal coordinate $x$ in the propagation direction. The criteria for variations in $c_p$ and $c_s$ that can be included under the mild non-uniformity approximation are discussed in supplementary material § 1.1.

Substitution of the expansions (2.12) and (2.14) into the interface (boundary) conditions of (2.5) (true for all $x$ and $t$) leads to the relations (cf. Hasselmann Reference Hasselmann1963)

(2.18)\begin{gather} \left.\begin{gathered} -\lambda(k^2 + \alpha^2)({\rm d}A_1\cos\alpha h + {\rm d}A_2\sin\alpha h) = 0,\quad {\rm at}\ y = h, \\ \alpha\,{\rm d}A_2 +{\rm i}\bar{\alpha}\,{\rm d}A' + {\rm i} k\,{\rm d}B' = 0,\quad {\rm at}\ y = 0, \\ -\lambda (k^2 + {\alpha}^2)\,{\rm d}A_1 + [\lambda' (k^2 + \overline{\alpha}^2) + 2\mu'\overline{\alpha}^2]\,{\rm d}A'_1 + 2\mu' k \bar{\beta}\,{\rm d}B_1'= 0,\quad {\rm at}\ y = 0, \\ -2{\rm i} k \bar{\alpha}\,{\rm d}A' + (k^2 - \bar{\beta}^2)\,{\rm d}B'_1 = 0,\quad {\rm at}\ y = 0. \end{gathered}\right\} \end{gather}

In other words, a vector ${\boldsymbol {\textrm {d}A}}$ comprised of the four Fourier coefficients in (2.18) would represent a free wave field when

(2.19)\begin{equation} \boldsymbol{\mathsf{D}}\,{\boldsymbol{{\rm d}A}} = \boldsymbol{0}, \end{equation}

where the matrix $\boldsymbol{\mathsf{D}}$ contains terms in $k$, $\omega$, $h$ and the Lamé constants $\lambda$, $\lambda '$ and $\mu '$, and is given by

(2.20)\begin{equation} \boldsymbol{\mathsf{D}} = \begin{bmatrix} -\lambda \gamma^2\cos \alpha h & -\lambda \gamma^2\sin \alpha h & 0 & 0\\ 0 & \alpha & {\rm i}\bar{\alpha} & {\rm i}k\\ -\lambda \gamma^2 & 0 & \lambda'(k^2 +\overline{\alpha}^2)+ 2\mu'{\overline \alpha}^2 & 2\mu' k \bar{\beta}\\ 0 & 0 & -2k\bar{\alpha} & k^2-\bar{\beta}^2 \end{bmatrix}, \end{equation}

where

(2.21)\begin{equation} \gamma^2 = {\left(\frac{\omega}{c_1}\right)}^2 = k^2 + \alpha^2. \end{equation}

The condition for non-trivial solutions for the potentials,

(2.22)\begin{equation} \Vert \boldsymbol{\mathsf{D}}\Vert = 0, \end{equation}

leads to algebraic relations that can be used to obtain dispersion surfaces on which $\omega$ and ${\boldsymbol k}$ for a particular propagation direction lie. Equation (2.22) has multiple roots, with each surface $(\omega _n, {\boldsymbol k}_n)$ defining the $n$th natural mode for the two-layer system under the present interface conditions.

Some special cases arise when $k = 0$, such as would occur when the two interacting surface waves producing the acoustic-gravity wave have equal wavelengths and opposite directions (see § 5) – for instance, when surface-wave conditions are influenced by a highly reflective shoreline. The acoustic-gravity waves then travel only vertically and generate a disturbance into the seafloor. In this case, the matrix $\boldsymbol{\mathsf{D}}$ reduces to $\boldsymbol{\mathsf{D}}_v$, where

(2.23)\begin{equation} \boldsymbol{\mathsf{D}}_v = \begin{bmatrix} -\lambda \alpha^2\cos\alpha h & -\lambda \alpha^2\sin\alpha h & 0 & 0\\ 0 & \alpha & {\rm i}\bar{\alpha} & 0\\ -\lambda \alpha^2 & 0 & \lambda' \overline{\alpha}^2+2\mu'{\overline \alpha}^2 & 0\\ 0 & 0 & 0 & -\bar{\beta}^2 \end{bmatrix},\end{equation}

and $\textrm {d}B_1' = 0$ would mean that there is no distortional particle oscillation component in this case. On the other hand, if the $3\times 3$ determinant multiplying the term $-{\overline \beta }^2$ is zero, that is,

(2.24)\begin{equation} \begin{vmatrix} -\lambda\alpha^2\cos\alpha h & -\lambda\alpha^2\sin\alpha h & 0\\ 0 & \alpha & {\rm i}\bar{\alpha} \\ -\lambda \alpha^2 & 0 & \lambda'\overline{\alpha}^2+2\mu'\overline{\alpha}^2 \end{vmatrix} = 0, \end{equation}

then recognizing that $\lambda '+ 2\mu ' = C'_{11}$ and $\lambda = C_{11}$ represent the elastic moduli for the seafloor and water, respectively (cf. Graff Reference Graff1991, chapter 6), (2.24) would give a relation between $\alpha$, $h$ and $\bar {\alpha }$ that includes the elastic moduli for water and the solid:

(2.25)\begin{equation} \alpha \tan \alpha h ={-}{\rm i}\,\frac{C'_{11}}{C_{11}}\,\bar{\alpha}. \end{equation}

With appropriate changes in notation and substitutions, (2.25) is seen to match (5.81) in Jensen et al. (Reference Jensen, Kuperman, Porter and Schmidt2011) for the vertical part of their solution for the compression wave (specifically, with $\bar {\alpha } \equiv k_{z,b}$, $\alpha \equiv k_z$ and $C_{11}'/C_{11} = \rho ' \alpha ^2/\rho \overline {\alpha }^2$), which admits dissipative leaky modes, thereby improving near-field accuracy of the underwater wave field (Jensen et al. Reference Jensen, Kuperman, Porter and Schmidt2011). The resulting improvement in near-field accuracy could be of consequence to ambient noise and energy conversion directly under regions where shoreward near-shore surface gravity waves interact with their reflections.

In general, the bulk waves and leaky modes may exist alongside seafloor surface-wave modes; here, the seafloor surface modes are studied further, given their closer relationship with microseisms. For seafloor surface waves, potential amplitudes decrease exponentially with depth below the seafloor surface (i.e. $y<0$). Using the substitutions

(2.26a,b)\begin{equation} \alpha' ={-}{\rm i}\bar{\alpha}\quad {\rm and}\quad \beta' ={-}{\rm i}\bar{\beta} \end{equation}

in (2.15), we have

(2.27)\begin{equation} \phi'(x,y; t) = \iint_{k, \omega} {\rm d}A'_1\,{\rm e}^{\alpha' y} \exp({{\rm i}(kx - \omega t)}) \end{equation}

and

(2.28)\begin{equation} H'_z(x,y; t) = \iint_{k, \omega} {\rm d}B'_1\,{\rm e}^{\beta' y}\exp({{\rm i}(kx - \omega t)}). \end{equation}

With these substitutions, the matrix in (2.20) becomes

(2.29)\begin{equation} \boldsymbol{\mathsf{D}} = \begin{bmatrix} -\lambda \gamma^2\cos \alpha h & -\lambda \gamma^2\sin \alpha h & 0 & 0\\ 0 & \alpha & -\alpha' & {\rm i}k\\ -\lambda \gamma^2 & 0 & \lambda'(k^2 -{\alpha'}^2)-2\mu'{\alpha'}^2 & 2{\rm i}\mu' k \beta'\\ 0 & 0 & -2{\rm i}k\alpha' & k^2+{\beta'}^2 \end{bmatrix}. \end{equation}

The dispersion relations resulting from (2.22) now can be expressed in terms of $(\omega, {\boldsymbol k})$ with real-valued, positive $\alpha '$ and $\beta '$ related to $\omega$ and $k$ according to

(2.30)\begin{equation} \left.\begin{gathered} \alpha' = \sqrt{(k^2 - \omega^2/c^2_p)}, \\ \beta' = \sqrt{(k^2 - \omega^2/c^2_s)}. \end{gathered}\right\} \end{equation}

Returning to the case $k = 0$, now

(2.31)\begin{equation} \alpha \tan \alpha h = \frac{\lambda' + 2\mu'}{\lambda}\,\alpha'. \end{equation}

It is noted that for $k = 0$, $\alpha '$ is still real-valued and positive for $\alpha h > 0$, so the vertical oscillations within the seafloor still decay to zero exponentially. The vertical pressure oscillations on the seafloor in this case can still generate Rayleigh-type waves with the Earth acting as a single semi-infinite elastic medium that is excited over a region by vertical pressure oscillations within a heavy compressible medium (cf. Stoneley Reference Stoneley1926).

If the seafloor is assumed to be rigid, then $\lambda ',\mu ' \rightarrow \infty$ in (2.31), which would imply that for a finite $\alpha$,

(2.32)\begin{equation} \tan\alpha h \rightarrow \infty,\quad {\rm or}\quad \cos\alpha h = 0. \end{equation}

Equation (2.32) leads to the vertical resonance modes given by (cf. Longuet-Higgins Reference Longuet-Higgins1950)

(2.33)\begin{equation} \cos \alpha_n h = 0\quad \Rightarrow\quad \frac{\omega_n h}{c_1} = (2n+1)\,\frac{\rm \pi}{2},\quad n = 0, 1, 2, \ldots. \end{equation}

The acoustic-gravity wave frequencies $\omega _n$ are then

(2.34)\begin{equation} \omega_n = \left(n+\frac{1}{2}\right)\frac{{\rm \pi} c}{h},\quad n = 0, 1, 2, \ldots. \end{equation}

If the acoustic-gravity waves are created by a second-order interaction of like-frequency surface waves from opposite directions – so that with $\sigma$ denoting the angular frequency of surface waves, $\omega _n = 2 \sigma _n$ – then

(2.35)\begin{equation} \sigma_n = \left(n+\frac{1}{2}\right)\frac{{\rm \pi} c}{2h},\quad n = 0, 1, 2, \ldots. \end{equation}

This conclusion is in agreement with the result derived by Longuet-Higgins (Reference Longuet-Higgins1950) for like-frequency waves from opposite directions over a rigid seafloor.

Returning to the original situation with a deformable seafloor, when $\beta ' = 0$ and $\textrm {d}B'_1 \neq 0$ in (2.23), the distortional potential $H_z$ would reduce to

(2.36)\begin{equation} H'_z = \int_{\omega} {\rm d}B'_1\exp({-{\rm i}\omega t}), \end{equation}

representing a non-propagating shear-mode oscillation within the seafloor.

Another special case arises when $\alpha = 0$. Referring back to (2.21), with $\alpha = 0$,

(2.37)\begin{equation} k = \frac{\omega}{c_1}, \end{equation}

and the compression waves travel only horizontally as plane waves, and span the entire water depth (cf. Kibblewhite & Wu Reference Kibblewhite and Wu1991). Since propagation then is in the horizontal direction only, the reflected component vanishes. In this case, the compression waves do not force the seafloor vertically. If they occur within canyons, then they may result in seiche modes in the water and seismic modes through the canyon walls, with the appropriate boundary conditions applied in the horizontal directions, though in the open ocean on a flat seafloor, if viscous effects are small, then they would not couple to seafloor waves.

When $k > \omega /c1$, $\alpha$ becomes imaginary, the downward component decays exponentially with depth, and the compression waves travel horizontally, but are confined to a region close to the water surface. Moreover, for $k \geq \omega /c_1$, free waves in water and in solid can exist independently of each other. Explanatory diagrams showing the different cases of interactions between surface wavenumbers and resulting pressure-wave wavenumbers are discussed by Kibblewhite & Wu (Reference Kibblewhite and Wu1991).

Within the wavenumber range $0 < k < \omega /c_1$, the compression wave paths are inclined, and their velocities have vertical and horizontal components. They undergo reflection at the seafloor, and some of their energy is transferred to the Scholte waves that they generate. However, the relations in (2.30) further restrict the wavenumber range for Scholte wave generation, since $\alpha '$ and $\beta '$ are real only if $\omega /c_p < k$ and $\omega /c_s < k$, respectively. Because $c_s \leq c_p$ typically, the admissible wavenumber range for the acoustic-gravity–Scholte wave system could be expressed as $\omega /c_s < k < \omega /c_1$, but more generally, it is needed that

(2.38)\begin{equation} \frac{\omega}{c_s}, \frac{\omega}{c_p} < k < \frac{\omega}{c_1}. \end{equation}

Over the $(\omega, k)$ range defined by (2.38), the acoustic-gravity and Scholte waves travel horizontally along the $x$ axis as a group, with a vertical variation as expressed in (2.14) (and/or (2.12)).

3. Propagation in mildly non-uniform media

In the treatment below, the dispersion relations for the inclined acoustic-gravity–Scholte wave coupled systems that result from (2.22) are cast in a form that is used next to study propagation in media with mild spatial non-uniformity (see Whitham (Reference Whitham1973, chapter 11) and supplementary material § 1.1). In particular, acoustic speed $c_1$ in water is taken to be a known function of the depth coordinate $y$, and the two propagation speeds $c_p$ and $c_s$ in the seafloor are assumed to change mildly with horizontal position ${\boldsymbol x}$. Water depth $h$ is still assumed to be uniform in this study. Free waves are considered first.

3.1. Free-wave propagation

From (2.22) and in terms of the pressure wave Fourier–Stieltjes coefficients,

(3.1a–c)\begin{equation} {{\rm d}}A_2 = f_1\,{\rm d}A_1,\quad {\rm d}A'_1 = f_2\,{\rm d}A_1\quad {\rm and}\quad {\rm d}B'_1 = f_3\,{\rm d}A_1, \end{equation}

where

(3.2)\begin{equation} \left.\begin{gathered} f_1 ={-}\cot \alpha h, \\ f_2 ={-}\frac{\alpha}{\alpha' - {\rm i} k C/D}\,f_1, \\ f_3 ={-}\frac{C}{D}\,f_2. \end{gathered}\right\} \end{equation}

It is noted that an alternative expression for $f_1$ can be stated:

(3.3)\begin{equation} f_1 = \left(\frac{-\lambda \gamma^2}{A +BC/D}\right) \div \left(\frac{-\alpha}{\alpha' + {\rm i}kC/D}\right). \end{equation}

Essentially, setting the two expressions for $f_1$ to be equal is equivalent to expanding the determinant $\Vert \boldsymbol{\mathsf{D}}\Vert = 0$. The quantities $A$, $B$, $C$ and $D$ are given by

(3.4)\begin{equation} \left.\begin{gathered} A = \lambda'(k^2 - {\alpha'}^2)-2\mu'{\alpha'}^2, \\ B = 2{\rm i}\mu'k\beta', \\ C ={-}2 {\rm i}k \alpha',\\ D = {\beta'}^2 + k^2. \end{gathered}\right\} \end{equation}

Now, $\textrm {d}A_1$ represents one of the Fourier coefficients of the vertically trapped acoustic-gravity wave field in the water column. It is recalled that the acoustic speed in water is related to $\rho$, the density of seawater, and $\lambda$, the Lamé constant (elastic modulus) of seawater. The seafloor Lamé constants $\lambda '$ and $\mu '$,along with the material density $\rho '$, give the two velocities $c_p$ and $c_s$. Summarizing,

(3.5a–c)\begin{equation} c^2_1 = \frac{\lambda}{\rho},\quad c^2_p = \frac{\lambda'+2\mu'}{\rho'}\quad {\rm and}\quad c^2_s = \frac{\mu'}{\rho'}, \end{equation}

where $c_1$ and $(c_p, c_s)$ are taken to be mildly varying functions of $y$ and $x$, respectively (see supplementary material § 1.1). The wavenumber is now included as a vector ${\boldsymbol k} = (k_1, k_2)$ in order to bring out the propagation direction $\theta _p$. For the given media (seawater and seafloor), $\alpha$, $\alpha '$ and $\beta '$ are determined for any ($\omega, k)$ on the dispersion surface for the given boundary conditions; $\alpha '$ and $\beta '$ represent solid particle motion that decays exponentially with depth beneath the seafloor. Specific wave groups can be tracked by following particular $(\omega, {\boldsymbol k})$ combinations on the dispersion surface. Individual wave groups travel at the group velocity $(C_{k_1}, C_{k_2})$, which remains constant for that group in uniform media, though it is a function of position in spatially non-uniform media. In spatially non-uniform media, $\omega$ and ${\boldsymbol k}$ for a particular group also change with position, whereas they would remain constant in uniform media. The dispersion surface in non-uniform media with mildly changing properties can be described using $\omega = W({\boldsymbol k}, {\boldsymbol x}, y)$ (Whitham Reference Whitham1973). Following the general discussion in (Whitham Reference Whitham1973, chapter 11), for propagation of the acoustic-gravity–Scholte wave system, if $\theta _{\pm } = ({\boldsymbol k}\boldsymbol {\cdot }{\boldsymbol x} \pm \alpha y- \omega t)$, then the consistency conditions resulting from $\theta _{tx} = \theta _{xt}$, $\theta _{x_1 x_2} = \theta _{x_2 x_1}$, $\theta _{ty} = \theta _{yt}$, and so on, yield

(3.6)\begin{equation} \left.\begin{gathered} \frac{\partial k_i}{\partial t} + \frac{\partial \omega}{\partial k_j}\, \frac{\partial k_i}{\partial x_j} \pm \frac{\partial \omega}{\partial \alpha}\, \frac{\partial k_i}{\partial y}+ \frac{\partial \omega}{\partial x_i} = 0,\quad i, j = 1, 2, \\ \frac{\partial \alpha}{\partial t} \pm \left(\frac{\partial \omega}{\partial k_j}\, \frac{\partial \alpha}{\partial x_j} + \frac{\partial \omega}{\partial \alpha}\, \frac{\partial \alpha}{\partial y} + \frac{\partial \omega}{\partial y}\right) = 0,\quad j = 1, 2, \\ \frac{\partial k_i}{\partial x_j} - \frac{\partial k_j}{\partial x_i} = 0,\quad \pm \frac{\partial \alpha}{\partial x_i} -\frac{\partial k_i}{\partial y} = 0,\quad i, j = 1, 2, i \neq j, \end{gathered}\right\} \end{equation}

with summation implied over the index $j$ in the first two conditions of (3.6). The $+$ and $-$ signs indicate downward and upward propagation, respectively. The group velocity components are defined as

(3.7a,b)\begin{equation} \frac{\partial \omega}{\partial k_i} = C_{k_i}, i = 1, 2,\quad {\rm and}\quad \frac{\partial \omega}{\partial \alpha} = C_{\alpha}. \end{equation}

Equation (3.6) then can be rewritten as

(3.8)\begin{equation} \left.\begin{gathered} \frac{\partial k_1}{\partial t} + C_{k_1}\, \frac{\partial k_1}{\partial x_1} + C_{k_2}\, \frac{\partial k_1}{\partial x_2} \pm C_\alpha\, \frac{\partial k_1}{\partial y}={-}\frac{\partial \omega}{\partial x_1}, \\ \frac{\partial k_2}{\partial t} + C_{k_1}\,\frac{\partial k_2}{\partial x_1} + C_{k_2}\, \frac{\partial k_2}{\partial x_2} \pm C_\alpha\,\frac{\partial k_2}{\partial y} ={-}\frac{\partial \omega}{\partial x_2}, \\ \frac{\partial \alpha}{\partial t} \pm \left(C_{k_1}\,\frac{\partial \alpha}{\partial x_1} + C_{k_2}\,\frac{\partial \alpha}{\partial x_2} + C_\alpha\,\frac{\partial \alpha}{\partial y}\right) ={\mp}\frac{\partial \omega}{\partial y}. \end{gathered}\right\} \end{equation}

The term on the right-hand side is non-zero in each equation because $\omega$ has an explicit dependence on ${\boldsymbol x}$ and $y$ in the two non-uniform media here. The characteristic curves for this system can be described as (Whitham Reference Whitham1973; see also Hasselmann Reference Hasselmann1963)

(3.9)\begin{equation} \left.\begin{gathered} \frac{{\rm d}{\boldsymbol x}}{{\rm d} t} = \frac{\partial W}{\partial {\boldsymbol k}}, \\ \frac{{\rm d}{\boldsymbol k}}{{\rm d} t} ={-}\frac{\partial W}{\partial{\boldsymbol x}}, \\ {\rm with}\quad\omega = W({\boldsymbol k},\alpha, {\boldsymbol x}, y). \end{gathered}\right\} \end{equation}

Because of the non-uniformity in the media, ${\boldsymbol k}$ and ${\omega }$ for a wave group vary with position as the group propagates, and the characteristics are not straight lines.

An ‘average variational principle’ of Whitham (Reference Whitham1973) could be used to study propagation in the ${\boldsymbol k}$ direction. For the random wave fields being studied this work, an average variational principle is proposed that utilizes the Fourier coefficients directly. Since the later steps lead to the evolution of wavenumber–frequency spectra, such an approach seems reasonable. Thus

(3.10)\begin{equation} J = \iiint_{t,x_1,x_2}\iint_{{\boldsymbol k},\omega}{\rm d}{\mathcal{L}} (\omega, k_1, k_2, {\rm d}A_1)\exp({{\rm i}{\boldsymbol k}\boldsymbol{\cdot} {\boldsymbol x} - \omega t})\,{\rm d} t\,{\rm d} x_1\,{\rm d} x_2, \end{equation}

where the integrals over each variable are from $-\infty$ to $\infty$. Here, $\textrm {d}{\mathcal {L}}$ is defined as

(3.11) \begin{equation} {{\rm d}}{\mathcal{L}} = G(\omega, {\boldsymbol k}, \alpha, \alpha', \beta', {\boldsymbol x}, y)\,{\rm d}A_1\,{\rm d}A^\ast_1, \end{equation}

and $\omega$ itself is $\omega = W({\boldsymbol k}, \alpha, \alpha ', \beta ', {\boldsymbol x}, y)$. Since the non-uniformities are mild in the sense of supplementary material § 1.1, an expression for $G$ is proposed as

(3.12)\begin{align} G &= [(\omega^2 - c^2_1(y)\,\alpha^2 - c^2_1(y)\,\vert{\boldsymbol k}\vert^2) + (\omega^2 - c^2_1(y)\,\alpha^2 - c^2_1(y)\,\vert{\boldsymbol k}\vert^2)\,F_1(\alpha)] \nonumber\\ &\quad + [(\omega^2 - c^2_p({\boldsymbol x})\,\alpha'^2 - c^2_p({\boldsymbol x})\, \vert{\boldsymbol k}\vert^2)\,F_2({\boldsymbol k}, \alpha,\alpha',\beta') \nonumber\\ &\quad +(\omega^2 - c^2_s({\boldsymbol x})\,\beta'^2 - c^2_s({\boldsymbol x})\, \vert{\boldsymbol k}\vert^2)\,F_3({\boldsymbol k},\alpha,\alpha',\beta')]. \end{align}

Here, to stress that $f_1$, $f_2$ and $f_3$ are not constants, the functions $F_1,F_2, F_3$ are introduced such that

(3.13)\begin{equation} \left.\begin{gathered} F_1(\alpha) = f_1 f^\ast_1, \\ F_2({\boldsymbol k},\alpha,\alpha',\beta') = f_2f^\ast_2, \\ F_3({\boldsymbol k},\alpha,\alpha',\beta') = f_3 f^\ast_3. \end{gathered}\right\} \end{equation}

The goal of the variational approach here is to determine a function $\textrm {d}A_1(\omega, {\boldsymbol k})$ that minimizes the action variable as the wave group propagates horizontally. For the proposed variational approach based on $\textrm {d}{\mathcal {L}}$, the necessary conditions for the minimum are

(3.14)\begin{equation} \left.\begin{gathered} \frac{\partial{({\rm d}{\mathcal{L}})}}{\partial ({\rm d}A_1)} = 0,\quad \frac{\partial{({\rm d}{\mathcal{L}})}}{\partial ({\rm d}A^\ast_1)} = 0, \\ \frac{\partial}{\partial t}\left(\frac{\partial({\rm d}{\mathcal{L}})}{\partial \omega}\right) - \frac{\partial}{\partial x_i}\left(\frac{\partial ({\rm d}{\mathcal{L}})}{\partial k_i}\right) - \frac{\partial}{\partial y}\left(\frac{\partial ({\rm d}{\mathcal{L}})}{\partial \alpha}\right) = 0, \end{gathered}\right\} \end{equation}

where the repeated index $i$ indicates a summation over $i$. In addition, the solutions need to satisfy the present consistency conditions (3.6). For a non-trivial $\textrm {d}A_1$, the first of the necessary conditions (3.14) just gives

(3.15)\begin{equation} G(\omega,{\boldsymbol k},\alpha,\alpha',\beta',{\boldsymbol x},y) = 0, \end{equation}

noting that $\omega = W({\boldsymbol k},\alpha,\alpha ',\beta ',x,y)$. Differentiation of (3.15) with respect to $k_i$ and $\alpha$ leads to

(3.16)\begin{equation} \left.\begin{gathered} G_{\omega}\,\frac{\partial \omega}{\partial k_i} + G_{ki} = 0\quad \Rightarrow\quad G_{k_i} ={-}C_{k_i} G_\omega, \\ G_\omega\,\frac{\partial \omega}{\partial \alpha} + G_\alpha = 0\quad \Rightarrow\quad G_\alpha ={-}C_\alpha G_\omega. \end{gathered}\right\} \end{equation}

The subscripts $\omega$ and $k_i$ here indicate differentiation with the respective quantities. Similarly, differentiation with respect to $x_i$ ($i = 1, 2$) and $y$ provides

(3.17)\begin{equation} \left.\begin{gathered} G_{\omega}\,\frac{\partial \omega}{\partial x_i} + G_{x_i} = 0\quad \Rightarrow\quad \frac{\partial \omega}{\partial x_i} ={-}\frac{G_{x_i}}{G_\omega},\quad i = 1, 2, \\ G_{\omega}\,\frac{\partial \omega}{\partial y} + G_y = 0\quad \Rightarrow\quad \frac{\partial \omega}{\partial y} ={-}\frac{G_\alpha}{G_\omega}. \end{gathered}\right\} \end{equation}

For the function $G$ in (3.12), $G_{x_i}$ and $G_y$ can be evaluated as

(3.18)\begin{equation} \left.\begin{gathered} G_{xi} = \frac{\partial G}{\partial c_p}\,\frac{\partial c_p}{\partial x_i} + \frac{\partial G}{\partial c_s}\,\frac{\partial c_s}{\partial x_i}, \\ G_y = \frac{\partial G}{\partial c_1}\,\frac{\partial c_1}{\partial y}. \end{gathered}\right\} \end{equation}

Sound speed $c_1$ in the ocean water column may in general vary with the latitude, longitude and depth below the surface. In this paper, only the dependence on the depth coordinate $y$ is accounted for.

Now,

(3.19)\begin{gather} \left.\begin{gathered} \frac{\partial ({\rm d}{\mathcal{L}})}{\partial \omega} = G_{\omega}\,{\rm d}A_1\,{\rm d}A^\ast_1 + G\,\frac{\partial}{\partial \omega}({\rm d}A_1\,{\rm d}A^\ast_1) = G_{\omega}\,{\rm d}A_1\,{\rm d}A^\ast_1, \\ \frac{\partial ({\rm d}{\mathcal{L}})}{\partial k_i} = G_{k_i}\,{\rm d}A_1\,{\rm d}A^\ast_1 + G\,\frac{\partial}{\partial k_i} ({\rm d}A_1\,{\rm d}A^\ast_1) = G_{k_i}\,{\rm d}A_1\,{\rm d}A^\ast_1={-}G_\omega (C_{k_i}\,{\rm d}A_1\,{\rm d}A^\ast_1), \end{gathered}\right\} \end{gather}

since $G = 0$, by the necessary condition (3.15). The second of the necessary conditions (3.14) now can be expressed as

(3.20)\begin{equation} \frac{\partial}{\partial t}(G_\omega\,{\rm d}A_1\,{\rm d}A^\ast_1) + \frac{\partial}{\partial x_i}(G_\omega C_{k_i}\,{\rm d}A_1\,{\rm d}A^\ast_1) + \frac{\partial}{\partial y}(G_\omega C_\alpha\,{\rm d}A_1\,{\rm d}A^\ast_1) = 0. \end{equation}

The non-uniformity of the media is taken into account in the following steps. The first term in (3.20) can be expanded as

(3.21)\begin{align} & \frac{\partial}{\partial t}(G_\omega\,{\rm d}A_1\,{\rm d}A^\ast_1) = \frac{\partial G_\omega}{\partial t}\,{\rm d}A_1\,{\rm d}A^\ast_1 + G_\omega \left(\frac{\partial}{\partial t}\,{\rm d}A_1\,{\rm d}A^\ast_1\right) \nonumber\\ &\quad =\frac{\partial G_\omega}{\partial k_i}\,\frac{\partial k_i}{\partial t}\,{\rm d}A_1\,{\rm d}A^\ast_1+ \frac{\partial G_\omega}{\partial \alpha}\,\frac{\partial \alpha}{\partial t}\,{\rm d}A_1\,{\rm d}A^\ast_1 + G_\omega \left(\frac{\partial}{\partial t}\,{\rm d}A_1\,{\rm d}A^\ast_1\right). \end{align}

The second term in (3.20) implies, for each $x_i$ in the summation,

(3.22)\begin{align} \frac{\partial}{\partial x_i}(G_\omega C_{k_i}\,{\rm d}A_1\,{\rm d}A^\ast_1) &= \frac{\partial G_\omega}{\partial k_i}\,\frac{\partial k_i}{\partial x_i}\,C_{ki}\,{\rm d}A_1\,{\rm d}A^\ast_1+ \frac{\partial G_\omega}{\partial \alpha}\,\frac{\partial \alpha}{\partial x_i}\,C_{k_i}\,{\rm d}A_1\,{\rm d}A^\ast_1 \nonumber\\ &\quad +\frac{\partial G_\omega}{\partial x_i}\,C_{k_i}\,{\rm d}A_1\,{\rm d}A^\ast_1 + G_\omega\,\frac{\partial}{\partial x_i}(C_{k_i}\,{\rm d}A_1\,{\rm d}A^\ast_1). \end{align}

The third term on the right-hand side of (3.22) arises because of the explicit dependence of $G_\omega$ on ${\boldsymbol x}$. The third term in (3.20) becomes

(3.23)\begin{align} \frac{\partial}{\partial y}(G_\omega C_\alpha\,{\rm d}A_1\,{\rm d}A^\ast_1) &= \frac{\partial G_\omega}{\partial k_i}\,\frac{\partial k_i}{\partial y}\,C_\alpha\,{\rm d}A_1\,{\rm d}A^\ast_1 + \frac{\partial G_\omega}{\partial \alpha}\,\frac{\partial \alpha}{\partial y}\,C_\alpha\,{\rm d}A_1\,{\rm d}A^\ast_1 \nonumber\\ &\quad +\frac{\partial G_\omega}{\partial y}\,C_\alpha\,{\rm d}A_1\,{\rm d}A^\ast{+} G_\omega\, \frac{\partial}{\partial y}(C_\alpha\,{\rm d}A_1\,{\rm d}A^\ast_1). \end{align}

Now, combining the first term in the second line of (3.21) with the first terms on the right-hand side of (3.22) and (3.23), and utilizing the consistency conditions (3.8), we have

(3.24)\begin{align} & \frac{\partial G_\omega}{\partial k_i} \left(\frac{\partial k_i}{\partial t}\,{\rm d}A_1\,{\rm d}A^\ast_1+ C_{k_j}\, \frac{\partial k_i}{\partial x_j}\,{\rm d}A_1\,{\rm d}A^\ast_1 + \frac{\partial k_i}{\partial y}\,{\rm d}A_1\,{\rm d}A^\ast_1\right) \nonumber\\ &\qquad +\frac{\partial G_\omega}{\partial \alpha} \left(C_{k_j}\,\frac{\partial \alpha}{\partial x_j}\,{\rm d}A_1\,{\rm d}A^\ast_1 + \frac{\partial \alpha}{\partial t}\,{\rm d}A_1\,{\rm d}A^\ast_1+ \frac{\partial \alpha}{\partial y}\,C_\alpha\,{\rm d}A_1\,{\rm d}A^\ast\right) \nonumber\\ &\quad ={-}\frac{\partial G_\omega}{\partial k_i}\,\frac{\partial \omega}{\partial x_i} - \frac{\partial G_\omega}{\partial \alpha}\,\frac{\partial \omega}{\partial y}. \end{align}

The final form of the second of the necessary conditions (3.20) is, for a non-zero $G_\omega$,

(3.25)\begin{align} & \frac{\partial}{\partial t}({\rm d}A_1\,{\rm d}A^\ast_1) + \frac{\partial}{\partial x_i}(C_{k_i}\,{\rm d}A_1\,{\rm d}A^\ast_1) + \frac{\partial}{\partial y}(C_\alpha\,{\rm d}A_1\,{\rm d}A^\ast_1) \nonumber\\ &\quad +\frac{1}{G_\omega}\left(-\frac{\partial G_\omega}{\partial k_i}\, \frac{\partial \omega}{\partial x_i}\,{\rm d}A_1\,{\rm d}A^\ast_1 + \frac{\partial G_\omega}{\partial x_i}\,C_{k_i}\,{\rm d}A_1\,{\rm d}A^\ast_1\right)\nonumber\\ &\quad +\frac{1}{G_\omega}\left(-\frac{\partial G_\omega}{\partial \alpha}\, \frac{\partial \omega}{\partial y}\,{\rm d}A_1\,{\rm d}A^\ast_1 + \frac{\partial G_\omega}{\partial y}\,C_\alpha\,{\rm d}A_1\,{\rm d}A^\ast_1\right) = 0. \end{align}

A wavenumber spectrum in terms of the propagating component ${\boldsymbol k}$ can be written as (Hasselmann Reference Hasselmann1963; Phillips Reference Phillips1977, chapter 3)

(3.26)\begin{equation} S_A({\boldsymbol k},t) = \frac{{\rm d}A_1\,{\rm d}A^\ast_1}{{\rm d}{\boldsymbol k}},\quad {\rm d}{\boldsymbol k} = {\rm d}k_1\,{\rm d}k_2. \end{equation}

Hence, dividing both sides of (3.25) by $\textrm {d}{\boldsymbol k}$, a relationship for the free-wave wavenumber spectrum for the acoustic-gravity waves propagating in the mildly non-uniform two-media system can be found:

(3.27)\begin{align} & \frac{\partial}{\partial t} S_A({\boldsymbol k},t) + \frac{\partial}{\partial x_i}(C_{k_i}\,S_A({\boldsymbol k},t)) + \frac{1}{G_\omega}\left(-\frac{\partial G_\omega}{\partial k_i}\, \frac{\partial \omega}{\partial x_i}\,S_A({\boldsymbol k},t) + \frac{\partial G_\omega}{\partial x_i}\,C_{k_i}\,S_A({\boldsymbol k},t)\right) \nonumber\\ &\quad + \frac{\partial}{\partial y}(C_\alpha\,S_A({\boldsymbol k},t)) + \frac{1}{G_\omega}\left(-\frac{\partial G_\omega}{\partial \alpha}\, \frac{\partial \omega}{\partial y}\,S_A({\boldsymbol k},t) + \frac{\partial G_\omega}{\partial y}\,C_\alpha\,S_A({\boldsymbol k},t)\right)= 0. \end{align}

Equation (3.27) describes how the wavenumber spectrum describing the potential $\phi$ evolves for free waves. Here, $S_A$ is a measure of energy over a small ball around $(\omega,{\boldsymbol k})$, and although $S_A$ remains constant along a group line in the absence of dissipation, the energy spreads horizontally over time. Equation (3.27) thus also represents the energy conservation (balance) relation for the wave system. The terms within the large parentheses represent the effect of the non-uniformities in the media, through $G_\omega$ and its derivatives with respect to $x_i$ and $y$, as well as through the $x_i$ and $y$ derivatives of the frequency $\omega$. As the wave group $(\omega,{\boldsymbol k})$ propagates in the ${\boldsymbol k}$ direction, both $\omega$ and ${\boldsymbol k}$ change along the group line. The group velocity components $C_{k_i}$ also change along the group line. This contrasts with propagation in uniform media, for which $(\omega,{\boldsymbol k})$ and $C_{k_i}$ for a group remain unchanged along the group line. The changes in group velocity signify continuous refraction of the acoustic-gravity–Scholte wave system as it propagates through the two media.

Equation (3.27) can be developed further by arguing that the partial derivative with respect to time in the first term, together with the partial derivatives with respect to $x_i$ and $y$, represents the material derivative of $S_A(\omega,{\boldsymbol k})$ along the group trajectory. In addition, $G_{\omega k_i} = G_{\omega \omega }({\partial \omega }{/\partial k_i})$, $\partial \omega /\partial k_i = C_{ki}$, and so on. Hence

(3.28) \begin{align} & \frac{{\rm d} S_A({\boldsymbol k},t)}{{\rm d} t} + \frac{1}{G_{\omega}} \left[{-}G_{\omega \omega}\,\frac{\partial \omega}{\partial x_i} + \frac{\partial G_\omega}{\partial x_i}\right]C_{ki}\,S_A({\boldsymbol k},t) \!+\! \frac{1}{G_{\omega}}\left[{-}G_{\omega \omega}\,\frac{\partial \omega}{\partial y} + \frac{\partial G_\omega}{\partial y}\right]C_\alpha\,S_A({\boldsymbol k},t) \nonumber\\ &\quad + S_A({\boldsymbol k},t)\left(\frac{\partial C_{ki}}{\partial x_i} + \frac{\partial C_{\alpha}}{\partial y}\right) = 0. \end{align}

Now, $G_\omega$, $G_{ki}$, $G_{\alpha }$, $C_{ki}$, $C_{\alpha }$, $G_{\omega \omega }$, $G_{xi}$, $G_y$, $\partial \omega /\partial x_i$ and $\partial \omega /\partial y$ can be expressed analytically (and evaluated in computer code) from (3.12), (3.16), (3.17), (3.18) and the assumed acoustic speed variations in (6.1). Further, for mildly non-uniform media (Whitham Reference Whitham1973),

(3.29)\begin{gather} C_{ki} = \frac{x_i}{t}, \end{gather}

(3.30)\begin{gather}\frac{\partial C_{ki}}{\partial x_i} = \frac{1}{t} + \frac{{\bar\partial} C_{ki}}{{\bar\partial} x_i}, \end{gather}

where $\bar \partial$ defines rate of change in $C_{ki}$ along the $x_i$ direction due to non-uniformity in the media. Similarly,

(3.31)\begin{equation} \frac{\partial C_\alpha}{\partial y} = \frac{1}{t} + \frac{{\bar\partial} C_\alpha}{{\bar\partial} y}, \end{equation}

where the overbar defines rate of change in $C_\alpha$ due to vertical non-uniformity. It is evident that for small $t$, the first term in (3.30) and (3.31) dominates other terms associated with non-uniformities, and $S_A$ varies as in uniform media. The effect of non-uniformities becomes more obvious as $t$ increases. Whereas in the horizontal plane the acoustic-gravity–Scholte waves can propagate to $x_i \rightarrow \infty$, the acoustic-gravity waves are confined within $0\leq y\leq h$ in the $y$ direction, undergoing a number of reflections off the two boundaries as the system propagates horizontally. Since some of the refraction effects during downward travel maybe undone during upward travel (the media properties are assumed to remain constant over time), an averaged quantity is used in the present calculations to represent the net effect of $C_\alpha$ dependence on $y$.

The Fourier coefficients $\textrm {d}A_1$ can be recovered from $S_A$ using the relation $\textrm {d}A_1 = \sqrt {S_A(\omega,{\boldsymbol k})\,\textrm {d}{\boldsymbol k}\,\textrm {d}\omega }$. The Fourier coefficients for acoustic-gravity pressure are related to $\textrm {d}A_1$ as

(3.32)\begin{equation} {\rm d}P_1 = \rho \omega^2\,{\rm d}A_1,\quad {\rm or}\quad {\rm d}P_1 = \lambda \gamma^2\,{\rm d}A_1, \end{equation}

where $\gamma = \omega /c_1 = \sqrt {k^2 + \alpha ^2}$. From (3.32), the wavenumber spectrum $S_P({\boldsymbol k})$ for the pressure variation in the water column can be found as

(3.33)\begin{equation} S_P({\boldsymbol k},t) = \frac{{\rm d}P_1\,{\rm d}P^\ast_1}{{\rm d}{\boldsymbol k}} = (\lambda \gamma^2)^2\,S_A({\boldsymbol k},t). \end{equation}

The wavenumber spectra $S_P$ can be used as estimates of low-frequency ambient noise spectra. The wavenumber spectrum (convertible to frequency spectrum) for the power density is then

(3.34)\begin{equation} S_W({\boldsymbol k},t) = \frac{S_P({\boldsymbol k},t)}{\rho c_1}. \end{equation}

The areal density of power $\varPi (\omega,{\boldsymbol k})$ for the $(\omega,{\boldsymbol k})$ group on the dispersion surface can be evaluated as

(3.35)\begin{equation} \varPi({\boldsymbol k},t) = \sqrt{S_W({\boldsymbol k},t)\,{\rm d}{\boldsymbol k}}. \end{equation}

Here, $\varPi$ is in units such as watts per square metre. Although the wavenumber–frequency spectra are used directly in calculations in this paper, the magnitude of the Fourier coefficient $\textrm {d}P$ gives an idea of second-order pressure amplitudes:

(3.36)\begin{equation} \vert {\rm d}P \vert = \sqrt{S_P({\boldsymbol k},t)\,{\rm d}{\boldsymbol k}},\quad {\rm d}{\boldsymbol k} = {\rm d}k_1\,{\rm d}k_2. \end{equation}

Next, the wavenumber–frequency spectrum for the Scholte waves forming part of this group can be found as,

(3.37)\begin{equation} \left.\begin{gathered} S_R({\boldsymbol k},t) = [\gamma^2_p\,F_2({\boldsymbol k},\alpha,\alpha',\beta) + \gamma^2_s\,F_3({\boldsymbol k},\alpha,\alpha',\beta')]\,\frac{{\rm d}A_1\,{\rm d}A^{{\ast}}_1}{{\rm d}{\boldsymbol k}} \\ = \omega^2\left[\frac{F_2({\boldsymbol k},\alpha,\alpha',\beta')}{c^2_p} + \frac{F_3({\boldsymbol k},\alpha,\alpha',\beta')}{c^2_s}\right]S_A({\boldsymbol k},t), \\ \vert {\rm d}R \vert = \sqrt{S_R({\boldsymbol k},t)\,{\rm d}{\boldsymbol k}}, \end{gathered}\right\} \end{equation}

where $\textrm {d}R$ represents the Fourier coefficient for the Scholte wave displacement amplitudes. Here, the relations (3.1a–c), (3.2) and (3.13) have been used, with $\gamma _p = \sqrt {k^2 - \alpha '^2}$ and $\gamma _s = \sqrt {k^2 - \beta '^2}$. Also, $S_R(\omega,{\boldsymbol k})$ is the wavenumber–frequency spectrum for Scholte wave amplitudes, which can be evaluated using the solution for $S_A(\omega,{\boldsymbol k})$; and $S_R$ is derived from Fourier coefficients $\textrm {d}A'_1$ and $\textrm {d}B'_1$, which satisfy both the kinematic and kinetic boundary conditions on the seafloor, in which $\textrm {d}A_1$ and $\textrm {d}A_2$ also play a part (see (2.18)). It is noted that a single formulation leading to the solution to (3.27) here provides the wavenumber spectra for the pressure in the compression waves (also representing ambient acoustic noise) (3.33), the power density (3.35), and the wavenumber spectra for the seafloor Scholte waves (3.37), with the non-uniformities assumed to be mildly varying.

3.1.1. Alternative approach to deriving $S_A$ and $S_R$

The wavenumber spectra for acoustic-gravity and Scholte wave amplitudes can be derived following a shorter pathway that albeit does not elucidate the effects of non-uniformity of the media as clearly as the more detailed approach above. Equation (3.20) can be expanded by treating the term $G_\omega \,\textrm {d}A_1\,\textrm {d}A^\ast _1$ as a combined entity, so that

(3.38)\begin{align} & \frac{\partial}{\partial t}(G_\omega\,{\rm d}A_1\,{\rm d}A^\ast_1) + C_{k_i}\,\frac{\partial}{\partial x_i}(G_\omega\,{\rm d}A_1\,{\rm d}A^\ast_1) + C_\alpha\,\frac{\partial}{\partial y}(G_\omega\,{\rm d}A_1\,{\rm d}A^\ast_1) \nonumber\\ &\quad +\left(\frac{\partial C_{ki}}{\partial x_i} + \frac{\partial C_\alpha}{\partial y}\right)G_\omega\,{\rm d}A_1\,{\rm d}A^\ast_1 = 0. \end{align}

The repeated index $i$ implies summation. Here, $G_\omega$ can be derived as

(3.39)\begin{equation} G_\omega = 2\omega[1 + F_1(\alpha) + F_2(\alpha, k, \alpha', \beta') + F_3(\alpha, k, \alpha', \beta')],\quad k = \vert{\boldsymbol k}\vert.\end{equation}

Hence it is seen that a multiplication by $G_\omega$ allows the boundary conditions to be included directly in the following steps. Dividing (3.38) through by $\textrm {d}{\boldsymbol k}$, and introducing a wavenumber spectrum $S_{BC}({\boldsymbol k},t)$ such that

(3.40)\begin{equation} S_{BC}({\boldsymbol k},t) = \frac{G_\omega\,{\rm d}A_1\,{\rm d}A^\ast_1}{{\rm d}{\boldsymbol k}},\end{equation}

it is found that

(3.41)\begin{align} \frac{\partial}{\partial t}S_{BC}({\boldsymbol k},t) + C_{k_i}\, \frac{\partial}{\partial x_i} S_{BC}({\boldsymbol k},t) + C_{\alpha}\, \frac{\partial}{\partial y} S_{BC}({\boldsymbol k},t) + S_{BC}({\boldsymbol k},t) \left(\frac{\partial C_{ki}}{\partial x_i} + \frac{\partial C_\alpha}{\partial y}\right) = 0. \end{align}

The first three terms represent propagation of a ‘modified’ energy term and can be combined into a single material derivative as

(3.42)\begin{equation} \frac{{\rm d}}{{\rm d} t}S_{BC}({\boldsymbol k},t) + S_{BC}({\boldsymbol k},t) \left(\frac{\partial C_{ki}}{\partial x_i} + \frac{\partial C_\alpha}{\partial y}\right) = 0.\end{equation}

Partial derivatives of $C_{ki}$ and $C_\alpha$ are found as expressed in (3.30) and (3.31). The other effects of non-uniformity are embedded in the multiplication by $G_\omega (\omega, {\boldsymbol k}, \alpha, \alpha ', \beta ')$. For particular wave groups, $S_A({\boldsymbol k},t)$ can be recovered from $S_{BC}({\boldsymbol k},t)$ using

(3.43)\begin{equation} S_A({\boldsymbol k},t) = \frac{S_{BC}({\boldsymbol k},t)}{G_\omega}. \end{equation}

Other spectra such as $S_P$ and $S_R$, and quantities such as $\varPi$, can be found as outlined in (3.34), (3.37) and (3.35), and so on. The alternative approach could serve as a ‘check’ in calculations since there are fewer numerical operations involved in it.

3.2. Uniform media

The treatment above simplifies considerably when the two media are assumed to be uniform. In particular, the second and third terms of (3.27), the second terms of (3.30) and (3.31), and the third terms of (S21) and (S22) in supplementary material § 1.2, tend to zero. For propagation in the horizontal plane, (3.27) and (3.42) become

(3.44)\begin{equation} \left.\begin{gathered} \frac{{\rm d}}{{\rm d} t} S_A({\boldsymbol k},t) + \frac{S_A( {\boldsymbol k},t)}{t} = 0, \\ \frac{{\rm d}}{{\rm d} t}S_{BC}({\boldsymbol k},t) + \frac{S_{BC}( {\boldsymbol k},t)}{t} = 0. \end{gathered}\right\} \end{equation}

Here, in the absence of refraction, the propagation of acoustic-gravity waves can be accounted for in terms of cumulative distance travelled (multiples of $h$ due to successive reflections). Equations (3.44) are used in the present calculations to provide a comparison for the more general case of mildly non-uniform media.

4. Forced wave response

As before, the origin of the coordinate system is at the seafloor, the $x$ axis is along the propagation direction ${\boldsymbol k}$, and the $y$ axis points vertically upwards, water depth being assumed to be uniform. When the two-media system of §§ 2 and 3 is excited by a stationary random external pressure distribution at its top boundary $y = h$, the ensuing acoustic-gravity–Scholte wave field must be consistent with the boundary conditions at the top and bottom of the water column. These boundary conditions are effectively embedded in the relationships of (3.2), and forced wave wavenumber spectra for the compression and Scholte wave fields can be interrelated easily using the relations in (3.33) and (3.37). Impulsive forcing at the seafloor or stationary random pressure disturbances within or on the water column could be addressed by defining suitably the excitation term in the analysis that follows. However, other approaches, such as that of Renzi & Dias (Reference Renzi and Dias2014), could be better suited for direct time-domain solutions in the case of storm-generated hydroacoustic waves and surface waves, and the approach of Meza-Valle et al. (Reference Meza-Valle, Kadri and Ortega2023) may be preferred in the case when time-domain understanding of rogue waves is needed using their faster-propagating compression wave signatures.

Here, the displacement potential $\phi _e$ for the acoustic-gravity waves can be related to the pressure excitation at $y=h$ with the help of the linear forced-wave equation as,

(4.1)\begin{equation} \frac{\partial^2\phi_e}{\partial t^2} - c^2_1 \nabla^2\phi_e = \frac{p(x,h;t)}{\rho};\quad c^2_1 = \lambda/\rho. \end{equation}

Here, as before, $\rho$ is seawater density, and $\lambda$ is the Lamé constant for seawater. We consider $\rho$ and $\lambda$ (and hence $c_1$) to be mildly varying functions of the $y$ coordinate in the sense of supplementary material § 1.1, so that (4.1) still describes the response to first order. Here, $p(x,h;t)$ denotes the second-order pressure due to interacting surface waves or a sea–surface pressure disturbance due to an atmospheric phenomenon. Section 5 discusses an approach to determination of this pressure using the properties of interacting surface-wave properties, assuming the surface-wave fields to be stationary random. The pressure excitation can be expanded using the Fourier–Stieltjes integral as

(4.2)\begin{equation} p(x,h,t) = \iint_{k,\omega}{\rm d}P(\omega,k)\exp({{\rm i}({\boldsymbol k}\boldsymbol{\cdot} {\boldsymbol x}-\omega t)}), \end{equation}

where $\textrm {d}P(\omega,k)$ denotes the Fourier coefficient for the pressure excitation ($k = \vert {\boldsymbol k}\vert$). The forced wave potential $\phi _e(x,y;t)$ is next expanded as

(4.3)\begin{equation} \phi_e(x,y;t) = \iint_{k,\omega} {\rm d}B_e(y,\omega,k)\exp({{\rm i}({\boldsymbol k} \boldsymbol{\cdot}{\boldsymbol x}-\omega t)}). \end{equation}

An expansion for the vertical variation of $\textrm {d}B(y,\omega,k)$ in terms of the vertical eigenfunctions $\psi _n(y)$ that satisfy the boundary conditions (2.18) is proposed:

(4.4)\begin{equation} {\rm d}B_e = \sum_{n=1}^{\infty} {\rm d}A_{1en}(\cos\alpha_n y + f_{1n}\sin\alpha_n y). \end{equation}

As indicated in § 2, the first orthogonal mode here is of most interest. Therefore letting $\textrm {d}A_{1e}$ denote $\textrm {d}A_{1e1}$ for convenience,

(4.5)\begin{equation} {\rm d}B(y,\omega,k) \approx {\rm d}A_{1e}(\omega,k)(\cos\alpha y + f_1\sin\alpha y). \end{equation}

It follows that

(4.6)\begin{equation} \phi_e(x,y;t) = \iint_{k,\omega}({\rm d}A_{1e}\, (\cos\alpha y + f_1\sin\alpha y)\exp({{\rm i}({\boldsymbol k} \boldsymbol{\cdot}{\boldsymbol x}-\omega t)}).\end{equation}

Substitution of $\phi _{e}$ into (4.1) leads to

(4.7)\begin{align} & \iint_{k,\omega} -\omega^2\,{\rm d}A_{1e}\,(\cos\alpha y + f_1\sin \alpha y) \exp({{\rm i}({\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x}-\omega t)}) \nonumber\\ &\qquad + c^2_1\iint_{k,\omega} \gamma^2\,{\rm d}A_{1e}\,(\cos\alpha y + f_1\sin \alpha y) \exp({{\rm i}({\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x}-\omega t)}) \nonumber\\ &\quad =\frac{1}{\rho}\iint_{k,\omega}{\rm d}P(\omega,k) \exp({{\rm i}({\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x}-\omega t)}), \end{align}

for all $x$ and $t$ in the domain. Next, both sides are multiplied by $(\cos \alpha y + f_1\sin \alpha y)$. Integration from $0$ to $h$, and use of the orthogonality of the cosine and sine functions, helps to reduce the terms inside the double integral on the left-hand side to

(4.8)\begin{equation} -\omega^2(1+f_1)\,{\rm d}A_{1e}\,\frac{h}{2} + c^2_1(k^2 + \alpha^2)(1+f^2_1)\,{\rm d}A_{1e}\,\frac{h}{2}. \end{equation}

Since the exciting pressure distribution $\textrm {d}P(\omega, k)$ acts on the surface $y=h$ and is not dependent on $y$, the right-hand side can be expressed as

(4.9)\begin{align} & \frac{1}{\rho}\int_{0}^{h}\iint_{k,\omega}{\rm d}P(\omega,k)\,(\cos \alpha y + f_1\sin \alpha y) \exp({{\rm i}({\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x}-\omega t)})\,{{\rm d}y} \nonumber\\ &\quad =\frac{1}{\rho}\iint_{k,\omega}{\rm d}P(\omega,k) \left(\frac{1}{\alpha}\,[\sin \alpha h + f_1(1-\cos\alpha h)]\right) \exp({{\rm i}({\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x}-\omega t)}), \end{align}

where $k$ and $\alpha$ represent the first mode of the dispersion relation. Since

(4.10)\begin{equation} \frac{\lambda(k^2+\alpha^2)}{\rho} = \frac{\lambda}{\rho}\, \gamma^2_1 = c^2_1\gamma_1^2 = \omega^2_1, \end{equation}

the temporal dynamic condition describing the Fourier coefficients for the forced waves becomes

(4.11)\begin{equation} (-\omega^2 + \omega^2_1)(1+f^2_1)\,{\rm d}A_{1e}(\omega,k) = \frac{2\,{\rm d}P(\omega,k)\,[\sin \alpha h + f_1(1-\cos\alpha h)]}{\rho \alpha h}, \end{equation}

where $\gamma _1$ and $\omega _1$ are, respectively, the quantities $\sqrt {k^2+\alpha ^2}$ and $\omega$ belonging to the dispersion surface for mode 1. The Fourier coefficient $\textrm {d}A_{1e}$ describing the acoustic-gravity component of the forced wave system thus is

(4.12)\begin{equation} {\rm d}A_{1e}(\omega,k) = \frac{2 [\sin \alpha h + f_1(1-\cos\alpha h)]}{\rho \alpha h(1+f^2_1)}\, \frac{{\rm d}P(\omega,k)}{-\omega^2 + \omega^2_1}. \end{equation}

It is evident that close to resonance, even small excitations can lead to significant amplitudes in the water column and on the seafloor. Hence it is worthwhile to study the response near resonance in more detail. The total wave field at a frequency $\omega$ is a linear superposition of the free wave field at $\omega _1$ and the forced wave field at $\omega$:

(4.13)\begin{equation} {\rm d}A_{1}(t) = {\rm d}A_{1n}\exp({-{\rm i}\omega_1 t}) + \frac{D_1\,{\rm d}P\exp({-{\rm i}\omega t})}{-\omega^2 + \omega^2_1}, \end{equation}

where

(4.14)\begin{equation} D_1 = \frac{2 [\sin \alpha h + f_1(1-\cos\alpha h)]}{\rho \alpha h(1+f^2_1)}. \end{equation}

The system in (4.13) is reminiscent of an undamped single-degree-of-freedom mechanical oscillator, for which the oscillation amplitude increases as $t$ when the excitation frequency equals the natural frequency (e.g. see Rao Reference Rao2010, chapter 3). When excitation is random as in the present, acoustic-gravity wave case for which $\omega _1$ represents the first resonance mode, $\omega$ may approach $\rightarrow \omega _1$ for random lengths of time at random intervals, hence a steady increase in the amplitude of $\phi$ is not expected. Therefore, a relationship is sought instead between the wavenumber spectra for the excitation and for the response. In particular, the following representation is sought for the response,

(4.15)\begin{equation} \int_{k} {\rm d}A_1\,{\rm d}A^\ast_1\exp({{\rm i} {\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x}}), \end{equation}

while for the excitation side, the following form is used,

(4.16) \begin{equation} \iint_{k,\omega} {\frac{D^2_1\,{\rm d}P\,{\rm d}P^\ast}{(-\omega^2 + \omega^2_1)^2}}\, ({\rm e}^{-{\rm i}\omega t} - {{\rm e}}^{-{\rm i}\omega_1 t}) ({\rm e}^{{\rm i}\omega t } + {\rm e}^{{\rm i}\omega_1 t}) \exp({{\rm i} {\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x}})\,\delta(\omega-\omega_1), \end{equation}

as it leads to a real-valued expression for the acoustic-gravity spectrum and agrees with the result in Hasselmann (Reference Hasselmann1963). As $\omega \rightarrow \omega _1$, the integrand on the right-hand side without the delta function approaches a $0/0$ indeterminacy. On application of L'Hospital's rule via

(4.17)\begin{equation} \lim_{\omega\rightarrow \omega_1} \frac{\dfrac{{\rm d}}{{\rm d}\omega}[D^2_1\,{\rm d}P\,{\rm d}P^\ast\, ({\rm e}^{-{\rm i}\omega t} - {\rm e}^{-{\rm i}\omega_1 t}) ({\rm e}^{{\rm i}\omega t} + {\rm e}^{{\rm i}\omega_1 t})]}{\dfrac{{\rm d}}{{\rm d}\omega} [(-\omega^2 + \omega^2_1)^2]}, \end{equation}

the integrand becomes

(4.18)\begin{equation} \frac{-2 {\rm i} t D^2_1\,{\rm d}P\,{\rm d}P^\ast}{(-\omega^2 + \omega^2_1) ({-}4\omega)}\exp({{\rm i}{\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x}}). \end{equation}

This integrand has poles at $\omega = \pm \omega _1$ and at $\omega = 0$. Next, substituting (4.18) into the right-hand-side integral in (4.16) as

(4.19)\begin{equation} \iint_{k,\omega}\frac{-2 {\rm i} t D^2_1\,{\rm d}P\,{\rm d}P^\ast}{(-\omega^2 + \omega^2_1) ({-}4\omega)}\exp({{\rm i} {\boldsymbol k} \boldsymbol{\cdot}{\boldsymbol x}})\,\delta(\omega-\omega_1), \end{equation}

the limit as $\omega \rightarrow \omega _1$ can be evaluated using contour integration, where it can be seen that the integral along the circular path with radius approaching infinity converges to zero. Then using just the residue at $\omega = \omega _1$, we have

(4.20)\begin{equation} \left.\begin{gathered} R_1 = \lim_{\omega\rightarrow \omega_1} \frac{-{\rm i}t D^2_1\,{\rm d}P\,{\rm d}P^\ast(k,\omega)}{2\omega (\omega+\omega_1)} = \frac{-{\rm i}t D^2_1\,{\rm d}P\,{\rm d}P^\ast({\boldsymbol k},\omega_1)}{4\omega^2_1}, \\ 2{\rm \pi} iR_1= \frac{{\rm \pi} t D^2_1\,{\rm d}P\,{\rm d}P^\ast({\boldsymbol k},\omega_1)}{2\omega^2_1}. \end{gathered}\right\} \end{equation}

The excitation side is then seen to be

(4.21)\begin{equation} \int_{\boldsymbol k} \frac{{\rm \pi} tD^2_1}{2\omega^2_1\,{\rm d}\omega}\,{\rm d}P\, {\rm d}P^\ast(\boldsymbol{k,}\omega_1)\exp({{\rm i} {\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x}}) = \int_{\boldsymbol k} \frac{{\rm \pi} tD^2_1}{2\omega^2_1}\, \frac{{\rm d}P\,{\rm d}P^\ast}{{\rm d}{\boldsymbol k}}\exp({{\rm i} {\boldsymbol k} \boldsymbol{\cdot}{\boldsymbol x}}){\rm d}{\boldsymbol k}. \end{equation}

The left-hand side of (4.15) can be expressed as

(4.22)\begin{equation} \int_{\boldsymbol k} \frac{{\rm d}A_1\,{\rm d}A^\ast_1}{{\rm d}{\boldsymbol k}} \exp({{\rm i} ({\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x})}){\rm d}{\boldsymbol k}, \end{equation}

noting that

(4.23)\begin{equation} \frac{{\rm d}A_1\,{\rm d}A^\ast_1}{{\rm d}{\boldsymbol k}} = S_A({\boldsymbol k},t), \end{equation}

where $S_A$ is the wavenumber–frequency spectrum for the acoustic-gravity waves. It is interesting to observe how the wavenumber spectrum for the acoustic-gravity potential evolves in time. Because the wavenumber spectrum $S_E$ for the excitation pressure is related to $\textrm {d}P$ and $\textrm {d}P^\ast$ as

(4.24)\begin{equation} \frac{{\rm d}P\,{\rm d}P^\ast}{{\rm d}{{\boldsymbol k}}} = S_E(\omega_1, {\boldsymbol k}), \end{equation}

the acoustic-gravity wavenumber spectrum evolution in time can be related to the excitation side as

(4.25)\begin{equation} \int_{\boldsymbol k}S_A({\boldsymbol k},t) \exp({{\rm i} {\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x}})\,{\rm d}{\boldsymbol k} = \int_{{\boldsymbol k}} \frac{{\rm \pi} tD^2_1}{2\omega^2_1}\,S_E(\omega_1,{\boldsymbol k}) \exp({{\rm i} {\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol x}})\,{\rm d}{\boldsymbol k}. \end{equation}

The spectrum $S_A({\boldsymbol k},t)$ at resonance grows linearly with time as long as excitation $S_E(\omega _1,{\boldsymbol k})$ at the natural frequency $\omega _1$ is present. Thus

(4.26)\begin{equation} S_A({\boldsymbol k},t)|_{DS} = \frac{{\rm \pi} t D^2_1}{2\omega^2_1}\,S_E(\omega_1,{\boldsymbol k})|_{DS}. \end{equation}

The symbol $|_{DS}$ indicates that $k = \vert {\boldsymbol k}\vert$ lies on the dispersion surface $({\boldsymbol k},\omega )$. Equation (4.26) shows the rate at which energy is being added to the acoustic-gravity waves (and by extension, also Scholte waves) excited when a pressure spectrum excites a particular $({\boldsymbol k}, \omega _1)$ that lies on the dispersion surface. Multiple $({\boldsymbol k}, \omega _1)$ combinations that lie on the dispersion surface can be accounted for using linear superposition. Of the multiple dispersion surfaces that satisfy (2.19) and (2.20), only the fundamental dispersion surface (dispersion mode 1) is considered here for the depth range $\sim$1000–3500 m.

Independently, (4.26) can be generalized to represent the cumulative effect of multiple generation areas that particular wave groups on the dispersion surface may encounter along their propagation paths. For the simpler case of propagation in uniform media, we have

(4.27)\begin{equation} S_A({\boldsymbol k},t) = \sum_{n=1}^{N_G} \frac{{\rm \pi} (t_{ne} - t_{ns}) D^2_1}{2\omega^2_1}\, S_E(\omega_1,{\boldsymbol k}) + \int_{t_{ne}}^t \frac{S_A({\boldsymbol k},\tau)}{\tau}\,{\rm d}\tau. \end{equation}

Here, $t > t_{ne}$ denotes the present time, $t_{ns}$ is the time at which a group enters the $n$th generation area, $t_{ne}$ is the time when it leaves that area, $\tau$ is the time integration variable, and $N_G$ is the number of generation areas. Multiple generation areas along the propagation path could enable reinforcement of spectral amplitudes affected by (i) higher-order detuning effects (Kadri & Akylas Reference Kadri and Akylas2016), and (ii) changes in surface-wave spectral content due to nonlinearities (Komen & Hasselmann Reference Komen and Hasselmann1996). The duration under each generation area of projected length $L_n$ along the propagation direction can be estimated as

(4.28)\begin{equation} t _{ne}-t_{ns} = \frac{L_n}{C_k}, \end{equation}

where $C_k$ is the horizontal group velocity of the wave group. The compression wave pressure wavenumber spectra (representing ambient noise), seafloor power densities and microseism amplitude wavenumber spectra under forcing $S_E$ can be evaluated using (3.33), (3.35) and (3.37), respectively, with the forced response wavenumber spectra $S_A$ evaluated using (4.27).

5. Excitation by interacting surface waves

To derive the wavenumber–frequency spectrum $S_E$ that excites the acoustic-gravity–Scholte wave fields of § 4, the approach of Hasselmann (Reference Hasselmann1963) is adopted here. Although several alternatives estimating the surface-wave excited frequency spectra at the seafloor have been proposed since (e.g. Hughes Reference Hughes1976; Webb & Cox Reference Webb and Cox1986; Kibblewhite & Wu Reference Kibblewhite and Wu1991), the approach of Hasselmann (Reference Hasselmann1963) is the most general since it enables a full account of all wave periods and directions represented within available records of wave spectra. This is particularly important for multi-mode surface-wave spectra (e.g. Ochi Reference Ochi1998, chapter 7), and in wave conditions where swells and wind waves have overlapping periods and nearly opposing directions. Such wave fields may occur when waves from different parts of a storm's trajectory intersect. Figure 2 shows an overview of the action of interacting waves and acoustic-gravity and Scholte waves.

Figure 2.

Schematic overview of the dynamics, with the interacting surface waves on the free surface, the acoustic-gravity waves and the Scholte waves. The cross-bars on the arrows provide a notional indication of the wavelengths of the surface waves and the acoustic-gravity and Scholte waves.

Following Hasselmann (Reference Hasselmann1963), the wave profile $\zeta (x_1, x_2; t)$ is assumed to be a stationary random process. It may be expanded in terms of Fourier coefficients as

(5.1)\begin{equation} \zeta(x_1, x_2; t) = \iint_{{\boldsymbol \kappa}, \sigma} {\rm d} Z({\boldsymbol\kappa},\sigma) \exp({{\rm i}({\boldsymbol\kappa}\boldsymbol{\cdot}{\boldsymbol x} - \sigma t)}), \end{equation}

where ${\boldsymbol \kappa }$ and $\sigma$ represent, respectively, the wavenumber vector and the angular frequency of the surface-wave system. The velocity potential $\varGamma (x_1, x_2, y; t)$ can be represented as

(5.2)\begin{equation} \varGamma(x_1, x_2, y; t) = \iint_{{\boldsymbol \kappa},\sigma} {\rm d}E({\boldsymbol \kappa}, \sigma)\exp({{\rm i}({\boldsymbol \kappa} \boldsymbol{\cdot}{\boldsymbol x} - \sigma t)})\,{\rm e}^{\kappa z},\quad z = y-h, \kappa = \vert{\boldsymbol \kappa}\vert. \end{equation}

Here, $\zeta$ and $\varGamma$ are related by the kinematic and kinetic boundary conditions to be satisfied approximately at the undisturbed free surface (Hasselmann Reference Hasselmann1963). The second-order pressure due to the interaction of the wavenumbers ${\boldsymbol \kappa }$ at frequencies $\sigma$ contained within the wave field can be represented as

(5.3)\begin{equation} p({\boldsymbol x}; t) ={-}\rho ({\boldsymbol \nabla} \varGamma ({\boldsymbol x}, y; t))^2, \end{equation}

with

(5.4)\begin{equation} \frac{\partial \varGamma}{\partial t} ={-}g\zeta,\quad {\rm or}\quad -\frac{{\rm i}\sigma}{g}\,\varGamma = \zeta\quad \Rightarrow\quad \varGamma = \frac{ig\zeta}{\sigma}. \end{equation}

Hence representing

(5.5)\begin{align} p({\boldsymbol x}; t) ={-}\rho g^2\iint_{{\boldsymbol \kappa},\sigma} \iint_{{\boldsymbol \kappa}', \sigma'} {\rm d}Z\,{\rm d}Z'\,(\kappa \kappa' - {\boldsymbol \kappa}{\boldsymbol \kappa}')\exp[-{\rm i}(({\boldsymbol \kappa} + {\boldsymbol \kappa}')\boldsymbol{\cdot}{\boldsymbol x} - (\sigma + \sigma')t)]\,\frac{1}{\sigma \sigma'}, \end{align}

where the term $(\kappa \kappa ' - {\boldsymbol \kappa }{\boldsymbol \kappa })$ arises from the $(\boldsymbol \nabla \varGamma )^2$ operation (the scalar term from differentiation in the vertical direction, the vector term from differentiation in the horizontal plane), with the integration over ${\boldsymbol \kappa }$ and ${\boldsymbol \kappa }'$ allowing interactions over the entire range of wavenumbers to be considered, the wavenumber–frequency spectrum for the second-order pressure in (5.5) can be found by multiplying the right-hand side by the complex conjugate of itself. Next, knowing that the surface-wave wavenumber–frequency spectra can be expressed as

(5.6a,b)\begin{equation} S_\zeta({\boldsymbol \kappa},\sigma) = \frac{{\rm d}Z\,{\rm d}Z^\ast}{{\rm d}{\boldsymbol \kappa}\,{\rm d}\sigma}\quad {\rm and}\quad S_{\zeta'}({\boldsymbol \kappa}',\sigma') = \frac{{\rm d}Z'\,{\rm d}Z'^\ast}{{\rm d}{\boldsymbol \kappa}'\,{\rm d}\sigma'}, \end{equation}

the wavenumber–frequency spectrum for $p({\boldsymbol x}; t)$, which forms the excitation spectrum $S_E({\boldsymbol k}, \omega )$ of § 4, is found to be

(5.7) \begin{align} S_E({\boldsymbol k},\omega) &= \rho^2 g^4 \iint_{{\boldsymbol \kappa},\sigma} \iint_{{\boldsymbol \kappa}', \sigma'} S_\zeta({\boldsymbol \kappa},\sigma)\, S_{\zeta}'({\boldsymbol \kappa}',\sigma')\,(\kappa \kappa' - {\boldsymbol \kappa}{\boldsymbol \kappa}')^2 \nonumber\\ &\quad \times\frac{1}{(\sigma \sigma')^2}\,[\delta( {\boldsymbol k} - ({\boldsymbol \kappa} + {\boldsymbol \kappa}'))][\delta(\omega - (\sigma + \sigma'))]\,{\rm d}{\boldsymbol \kappa}\, {\rm d}{\boldsymbol \kappa}'\,{\rm d}\sigma \,{\rm d}\sigma'. \end{align}

Since only three waves ${\boldsymbol \kappa }$, ${\boldsymbol \kappa }'$ and ${\boldsymbol k}$ are involved in the interaction, the ${\boldsymbol k}$ wave must follow a dispersion relation different to that of the two surface waves ${\boldsymbol \kappa }$ and ${\boldsymbol \kappa }'$ for there to be meaningful energy transfer. This is possible when ${\boldsymbol k}$ is the horizontal wavenumber of an acoustic-gravity wave, which travels at a much greater phase velocity and follows a dispersion relation that allows such a three-wave interaction. In order to satisfy the acoustic-gravity dispersion relations, $k$ must therefore be much smaller than $\kappa$ and $\kappa '$, while $\omega$ is at least the same order of magnitude as $\sigma$. Therefore, it is mostly ${\boldsymbol \kappa }$ and ${\boldsymbol \kappa }'$ that are nearly opposed in directions that would dominate the integral in (5.7) (see schematic illustration in figure 3), in the process implying that the acoustic-gravity wave frequency would be nearly twice the surface-wave frequency. Although multiple combinations $({\boldsymbol \kappa }, {\boldsymbol \kappa }',\sigma, \sigma ')$ would be included within the integral in (5.7), an implicit further condition is that the relevant $({\boldsymbol \kappa },\sigma )$ combinations fall within the wavenumber and frequency ranges of practical surface-wave conditions. (See Phillips (Reference Phillips1960, Reference Phillips1977, chapter 3) for a discussion of the surface-wave dispersion relation precluding three-waveinteractions among three surface waves but allowing four-wave surface-wave interactions.)

Figure 3.

Schematic concept diagram illustrating the interaction between two surface waves with wavenumbers ${\boldsymbol \kappa }$ and ${\boldsymbol \kappa }'$ at frequencies $\sigma _1$ and $\sigma _2$. When two surface waves are nearly parallel, their interaction would lead to a horizontal wavenumber ${\boldsymbol k}$ that may be too large to lie on the dispersion surface for the acoustic-gravity–Scholte wave system (a), though with nearly opposing surface waves, the resulting horizontal wavenumber is small enough to lie on the dispersion surface (b), so resonant energy transfer can occur. The contribution of wave components like (b) will therefore dominate the integral in (5.7).

Here, as captured by the delta functions in (5.7), the wavenumbers and frequencies of the interacting three-wave system are related according to e.g. Hasselmann (Reference Hasselmann1966) and Ardhuin (Reference Ardhuin2018) as

(5.8a,b)\begin{equation} {\boldsymbol \kappa} + {\boldsymbol \kappa}' = {\boldsymbol k}\quad {\rm and}\quad \sigma + \sigma' = \omega. \end{equation}

It has been observed that the dominant surface-wave frequencies over long time periods are $\sim$0.1–0.11 Hz (or $\sigma,\sigma ' \sim 0.628$ rad s$^{-1}$), corresponding to swells of periods 9–10 s that travel long distances over the ocean surface after being generated by one or more storms (Webb Reference Webb1986). Therefore, it is plausible that the dominant acoustic-gravity–Scholte wave frequencies would be in the range 0.2–0.22 Hz (or $\omega \sim 1.26$ rad s$^{-1}$). Consequently, when two 10 s swells of wavenumbers $\kappa \sim 0.0162$ rad m$^{-1}$ interact in a depth of 2800 m for which the acoustic-gravity–Scholte wave dispersion relationship gives $k = 5.78 \times 10^{-4}$ rad m$^{-1}$ for a soft rock type seafloor, there is only a small range of swell directions, and the resulting acoustic-gravity–Scholte wave group direction for which the relationship ${\boldsymbol \kappa } + {\boldsymbol \kappa }' = {\boldsymbol k}$ will be satisfied, giving $\sigma + \sigma ' = \omega = 1.256$ rad s$^{-1}$. The $k$ value from the dispersion relation just implies resonant energy transfer. Hence when the surface-wave directional spectra are available, it is relatively straightforward to estimate, using a formal version of the diagram in figure 3, the direction in which bulk of the energy from the second-order interaction of the surface waves will travel as acoustic-gravity–Scholte waves.

These considerations are important for siting seafloor energy converter arrays in given geographic locations in the ocean, or when determining the angular ranges for the seafloor or coastal seismometers that would record particular storm activity at a particular time. As mentioned, a multi-directional wave field with overlapping wind-wave and swell conditions could be generated by a single storm if it travels over a curved trajectory.

Assuming that a multi-directional wave field is translating over the ocean surface at a velocity ${\boldsymbol U}$, which may be changing as the storm propagates (see, for instance, Marshall & Plumb Reference Marshall and Plumb2008, chapter 6), we have

(5.9)\begin{equation} \left(\frac{{\rm d}{\boldsymbol U}}{{\rm d} t}\right)_O = \left(\frac{{\rm d}{\boldsymbol U}}{{\rm d} t}\right)_o + \boldsymbol{\varOmega}\times {\boldsymbol U}, \end{equation}

where $O$ is the origin of an Earth-fixed reference frame, while $o$ is the origin of a reference frame travelling with the storm. Here, $\boldsymbol {\varOmega } = (0, 0, f_3)$, with $f_3 = f_e\cos L_t$ representing the Earth's angular velocity component at the latitude $L_t$ at which the storm-centred coordinate system origin $o$ is placed. The translation of the present acoustic-gravity–Scholte wave system is thus referenced to the storm-centred coordinate system $o- x_1 x_2 y$. The frequencies as expressed in $o- x_1 x_2 y$ therefore need to be adjusted according to the velocity ${\boldsymbol U}$ at which the surface-wave system is travelling (Phillips Reference Phillips1977, chapter 3) as expressed relative to the Earth-fixed reference frame at $O$.

For instance, if a single storm causes the two interacting surface waves of frequency $\sigma = 0.628$ rad s$^{-1}$ ($f = 0.1$ Hz), but currently the storm is moving at ${\boldsymbol U} = 10$ km h$^{-1}$ at an angle $\theta _u$, then the two interaction equations (5.8a,b) in the moving reference frame need to be amended as

(5.10)\begin{equation} \left.\begin{gathered} {\boldsymbol \kappa} + {\boldsymbol \kappa}' = {\boldsymbol k}, \\ \sigma + \sigma' = \omega', \\ \omega = \omega' + {\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol U}. \end{gathered}\right\} \end{equation}

The frequency recorded by a seafloor-fixed sensor will be $\omega$, while $\omega ' = 1.256$ rad s$^{-1}$. Since ${\boldsymbol k}\boldsymbol {\cdot }{\boldsymbol U} = 1.6\times 10^{-3}\cos (\theta _p-\theta _u)$, we have $\omega \in [1.2544, 1.2576]$ rad s$^{-1}$, as $(\theta _p-\theta _u) \in [0, {\rm \pi}]$. Greater differences may be observed for fast-moving storms.