%20-%20https://sketch.com%20--%3e%3ctitle%3e3ECBBF50-FD08-412A-943D-2CFC5E025A16%3c/title%3e%3cdesc%3eCreated%20with%20sketchtool.%3c/desc%3e%3cg%20id='Stylesheet'%20stroke='none'%20stroke-width='1'%20fill='none'%20fill-rule='evenodd'%3e%3cg%20transform='translate(-926.000000,%20-668.000000)'%20id='BUTTONS'%3e%3cg%20transform='translate(100.000000,%20597.000000)'%3e%3cg%20id='BTN-/-RIGHTS'%20transform='translate(816.000000,%2065.000000)'%3e%3cg%20id='ICON-/-Rights'%20transform='translate(10.000000,%206.000000)'%3e%3cpath%20d='M13.5752841,0%20L1.34943182,0%20C0.59375,0%200,0.586666667%200,1.33333333%20L0,22.6666667%20C0,23.4133333%200.620738636,24%201.34943182,24%20L17.6505682,24%20C18.40625,24%2019,23.3866667%2019,22.6666667%20L19,5.36%20L13.5752841,0%20Z%20M13.5752841,1.89333333%20L17.0838068,5.36%20L13.5752841,5.36%20L13.5752841,1.89333333%20Z%20M1.34943182,22.6666667%20L1.34943182,1.33333333%20L12.1988636,1.33333333%20L12.1988636,5.36%20C12.1988636,6.10666667%2012.7926136,6.69333333%2013.5482955,6.69333333%20L17.6235795,6.69333333%20L17.6505682,22.6666667%20L1.34943182,22.6666667%20Z'%20id='Shape'%20fill='%23595959'%20fill-rule='nonzero'%3e%3c/path%3e%3ccircle%20id='Oval'%20stroke='%23595959'%20stroke-width='1.5'%20fill='%23595959'%20cx='9.5'%20cy='13.5355339'%20r='5.5'%3e%3c/circle%3e%3cpolyline%20id='Path-2'%20stroke='%23FFFFFF'%20stroke-width='1.5'%20transform='translate(9.500000,%2013.035534)%20rotate(45.000000)%20translate(-9.500000,%20-13.035534)%20'%20points='8%2016.0355339%2011%2016.0355339%2011%2010.0355339'%3e%3c/polyline%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
1. Introduction
The evolution of the pressure field resulting from an initial disturbance in the ocean has been of significant interest to the scientific community for a long time, primarily due to the potential impact of underwater nuclear bomb testing, underwater volcanic eruptions, seismic events and the associated tsunami hazards. Due to these, there is widespread monitoring of ocean pressure; for example, the Comprehensive Nuclear-Test-Ban-Treaty Organization monitors pressure signals from hydrophone stations to detect possible nuclear tests, and similar devices are also used to detect tsunami waves. Accurate modelling of the resulting pressure signatures allows information about the initial event to be determined; for example, an attempt has been made to estimate the crash site location of Malaysian Airlines Flight 370 by solving an inverse problem (Kadri Reference Kadri2024). The accuracy of such inverse problems depends on the accuracy of the forward solution and analysis of the associated wave evolution problem. Obviously, the short–time behaviour of such events may be highly nonlinear and involve the formation of a cavity and subsequent expansion at the initial state. However, after some time (and at some distance), the pressure in the water column can be accurately modelled with a linearised water wave theory that accounts for the effect of dynamic compression (we use this term to distinguish the effect of the local compressibility from the static compression which is the accumulation of compression due to the weight of the overlying water column).
The inclusion of the compressibility of water is an essential part of our work. Although the compressibility of water is small, its inclusion changes the mathematical structure of the problem so that acoustic–gravity waves (AGW) arise. These are hydroacoustic waves influenced by gravity at low frequency, which are analogous to propagating modes in acoustic or electromagnetic waveguides. These modes capture the effect of the free-surface boundary condition and play a pivotal role in the accurate calculation of the pressure of the water column (Nosov et al Reference Nosov, Karpov, Kolesov, Sementsov, Matsumoto and Kaneda2018). It is worth mentioning that, in the limit of infinite sound speed in water, the gravity mode transforms into the propagating wave mode under realistic parameter choices, while the AGW modes become evanescent. Following the devastating aftermath of the
$2004$
Indian Ocean tsunami, a renewed interest in the application of AGW in tsunami detection mechanisms (through modelling of surface wave propagation) has grown through the works of Nosov (Reference Nosov, Karpov, Kolesov, Sementsov, Matsumoto and Kaneda1999), Nosov & Skachko (Reference Nosov and Skachko2001) and Nosov & Kolesov (Reference Nosov and Kolesov2007), which are based on the well-established fundamental mathematical theory of compressible ocean developed by Longuet-Higgins (Reference Longuet-Higgins1950), Miyoshi (Reference Miyoshi1954), Sells (Reference Sells1965) and Yamamoto (Reference Yamamoto1982). While incompressible water wave theory is most commonly used in practice, the importance of compressibility has been demonstrated by enhanced accuracy in estimating the observed tsunami wave profile (Nosov et al Reference Nosov, Karpov, Kolesov, Sementsov, Matsumoto and Kaneda2018). Following this development, Stiassnie (Reference Stiassnie2010) and Hendin & Stiassnie (Reference Hendin and Stiassnie2013) computed the far-field profiles of the gravity and AGW modes. Subsequently, there has been a surge in studies investigating different properties of these AGW modes, for example: the production of AGW by interacting gravity waves (Kadri & Stiassnie Reference Kadri and Stiassnie2013); triad resonances (Kadri Reference Kadri2016; Kadri & Akylas Reference Kadri and Akylas2016) – an interaction between two waves resulting in a third wave satisfying resonance conditions; their impact on different oceanic processes (e.g. deep-water transport (Kadri Reference Kadri2014), ice shelf breaking (Kadri & Stiassnie Reference Kadri and Stiassnie2012), wave dissipation due to sea ice (Chen, Gilbert & Guyenne Reference Chen, Gilbert and Guyenne2019)); and their response to different seafloor characteristics (Chierici, Pignagnoli & Embriaco Reference Chierici, Pignagnoli and Embriaco2010; Abdolali, Kirby & Bellotti Reference Abdolali, Kirby and Bellotti2015). Different mathematical techniques are also deployed to solve the associated initial boundary value problem (IBVP) (e.g. rigid but variable bathymetry using depth-integrated mild-slope equations (Sammarco et al Reference Sammarco, Cecioni, Bellotti and Abdolali2013), Fourier-series expansion (Das & Meylan Reference Das and Meylan2023b
) and eigenfunction expansion (Das & Meylan Reference Das and Meylan2023a
)), which are essential to address the static compression.
The static compression of the ocean is due to the effect of the slight compressibility of water, causing an increase in the density of the deep ocean. It is only significant for waves which penetrate the deep ocean, and even then, its effect is small. The initial account of static compression was given by Longuet-Higgins (Reference Longuet-Higgins1950). However, until the work of Abdolali & Kirby (Reference Abdolali and Kirby2017), showing the increase of tsunami phase speed at high ocean depth due to static compression, this quantity was rarely included in models due to its presumed insignificant effect on wave propagation. A few contemporary works retained static compression in the mathematical treatment of physical problems to study the behaviour of the gravity mode (Kadri & Stiassnie Reference Kadri and Stiassnie2013; Kadri Reference Kadri2015; Kadri Reference Kadri2016; Abdolali & Kirby Reference Abdolali and Kirby2017; Abdolali, Kadri & Kirby Reference Abdolali, Kadri and Kirby2019), most notably by Omira et al. (Reference Omira, Ramalho, Kim, González, Kadri, Miranda, Carrilho and Baptista2022) to match the Hunga Tonga–Hunga Ha‘apai volcanic eruption-induced tsunami data using the resonance between the acoustic waves and the surface gravity wave. Variations in water salinity will also affect the density of the ocean over its depth. As such, the vertical profile of sound speed can be complicated. While the common approach is to approximate the sound speed profile as constant, Michele & Renzi (Reference Michele and Renzi2020) used a perturbation technique on the variable-depth profile of the sound speed to derive an analytical expression correct up to third order that governs the evolution of the hydroacoustic waves. More recently, Andrade et al (Reference Andrade, Bustamante, Kadri and Stuhlmeier2025) used a dynamic approach to study triad resonances and energy exchange between different modes. Recent developments include a novel Lagrangian method to approach density stratification by linearising the compressible Euler equation (Dubois et al Reference Dubois, Imperiale, Mangeney, Bouchut and Sainte-Marie2023), who demonstrated that other models in the literature can be derived from their general proposed model. In the context of tsunami wave generation, Das & Meylan (Reference Das and Meylan2023a ) solved for the time-dependent evolution due to impulsive forcing from the sea floor, including the effect of static compression. The eigenfunction matching required a special inner product, which is generalised in the present work. This method was later extended to include more general geometries (e.g. a step (Das, Pethiyagoda & Meylan Reference Das, Pethiyagoda and Meylan2024), a trench (Pethiyagoda et al Reference Pethiyagoda, Das and Meylan2024b ), a generic ocean floor (Pethiyagoda et al Reference Pethiyagoda, Das and Meylan2025b )), three-dimensional solutions (Das & Meylan Reference Das and Meylan2024) and surface disturbance-induced wave propagation models (Pethiyagoda et al Reference Pethiyagoda, Das, Bonham and Meylan2024a , Reference Pethiyagoda, Das and Meylan2025a ). These studies augment the works where static compression is dropped (Renzi & Dias Reference Renzi and Dias2014) and ocean compressibility is ignored (Liu & Higuera Reference Liu and Higuera2022). We note here that the current work could be extended to include other effects that are non-standard beyond the static compression considered here. Two examples are the effect of the seabed elasticity (Eyov et al Reference Eyov, Klar, Kadri and Stiassnie2013; Abdolali et al Reference Abdolali, Kadri and Kirby2019; Kadri Reference Kadri2019; Williams & Kadri Reference Williams and Kadri2023) and more complex static compression using the vertical profile of the sound speed (Michele & Renzi Reference Michele and Renzi2020). We note that, for the range of parameters we investigate here, the effect of the static compression is small, even for the simplest formulation which we use.
Many energy-conserving linear wave problems can be recast as an abstract wave equation
$\partial _t^2 u+\mathcal{A}u=0$
, where
$\mathcal{A}$
is a positive self-adjoint operator in a Hilbert space
$(\mathcal{H},\langle \boldsymbol{\cdot },\boldsymbol{\cdot }\rangle )$
and
$u$
is the unknown variable. In one of the simplest cases, choosing
$\mathcal{A}=-c^2\partial _x^2$
with inner product
$L_2(\mathbb{R})$
gives the classical wave equation with wave speed
$c$
. More complicated forms of
$\mathcal{A}$
appear in various scenarios such as hydroelasticity (Meylan Reference Meylan2002). While the use of the spectral theorem to obtain the solution as an expansion over the eigenfunctions of
$\mathcal{A}$
was long known for the case when
$\mathcal{A}$
has a discrete spectrum (Fourier/Sturm Liouville theory) the case when
$\mathcal{A}$
has a continuous spectrum began being worked out for the Schrödinger equation with the advent of quantum mechanics (Povzner Reference Povzner1953; Ikebe Reference Ikebe1960). The theory, which has been extended to the classical wave equation (Wilcox Reference Wilcox1975) and to hydrodynamics for the case of incompressible water wave scattering (Hazard & Loret Reference Hazard and Loret2007; Hazard & Meylan Reference Hazard and Meylan2008; Meylan Reference Meylan2009), is known as a generalised eigenfunction expansion. The theoretical foundation is given by a Gelfand triple, also known as a rigged Hilbert space, which is a mathematical structure used to extend the usual Hilbert space framework to include more generalised functions (i.e. distributions). This allows the treatment of generalised eigenfunctions, i.e. functions
$\hat {u}$
satisfying
$\mathcal{A}\hat {u}=\omega ^2\hat {u}$
but which are not square-integrable (i.e. not belonging to the Hilbert space
$\mathcal{H}$
) where
$\omega$
is the angular frequency. These generalised eigenfunctions are still useful and meaningful in a generalised sense. In general, this theory has been developed in isolation from computation, despite the fact that the simplest example of this structure is the Fourier transform, and its computation by the discrete Fourier transform (often implemented using the fast Fourier transform) is one of the most ubiquitous algorithms.
It turns out that the theory of generalised eigenfunctions and spectral solutions for continuous spectra can be implemented numerically by applying quadrature to the integrals over the continuous spectrum. This idea was in the works of Hazard & Meylan (Reference Hazard and Meylan2008), Meylan (Reference Meylan2009), Meylan & Eatock Taylor (Reference Meylan and Eatock Taylor2009) and Peter & Meylan (Reference Peter and Meylan2010), but appears explicitly in Wilks et al (Reference Wilks, Meylan, Montiel and Wakes2024, Reference Wilks, Meylan, Montiel and Wakes2025). As far as we are aware, the present work is the first application of the theory to compressible water wave scattering. A key advantage of the numerical spectral technique is that the generalised eigenfunctions, which satisfy
$\mathcal{A}\hat {u}=\omega ^2\hat {u}$
, are precisely the time harmonic solutions for the angular frequency
$\omega$
. These single-frequency scattering solutions can be found numerically for a wide range of problems, e.g. variable bathymetry, which would be a relatively straightforward extension of the present work. We note that the ease of determining the spectrum stands in contrast with a more numerical method (e.g. Lanczos algorithm), in which the operator
$\mathcal{A}$
is approximated by a matrix and the spectrum of this matrix is sought. In this case, very many numerical problems arise, especially because the eigenvalues of the underlying operator are continuous. Further problems can arise from the fact that it may be difficult to preserve the self-adjoint property of the matrix if the inner product
$\langle \boldsymbol{\cdot},\boldsymbol{\cdot}\rangle$
is not straightforward to compute (as in the case with the current problem). Finally, we note that the use of these spectral methods is especially useful in problems where there is an initial condition rather than forcing at the boundary.
In this work, we consider the linear initial value problem for an underwater pressure source, including the effect of the free surface. The solution relies on the use of a special Hilbert space. We also extend this formulation to the more challenging case where we include the effect of static compression. The properties of the operator involved in this class of problems are not well studied in the literature. We provide a mathematical structure for the operators involved in a rigid and flat ocean floor via a special inner product. Additionally, the mathematical complexity entangled with static compression makes the problem a difficult one to solve in finite water depth. We aim to show the ease with which the present method can tackle such a complex problem, relaxing any assumptions on the relative water depth, thus making the approach a generalised one.
The manuscript is structured as follows. Section 2 presents the two-dimensional mathematical formulation of the physical problem where static compression is ignored, followed by the solution method that involves deriving a generic two-dimensional inner product in § 3. The time-domain simulations are detailed in § 4. The problem is much simpler if static compression is ignored, and for many scenarios, only this solution will be required, and this is the motivation for giving this in detail. In § 5, we extend the formulation to static compression and calculate the time-domain simulations of the flow field in this case. We note that the inclusion of static compression complicates the derivation, but the numerical solution takes a similar computational time to the solution without compression. The manuscript concludes with a brief summary in § 6.
2. Mathematical formulation without static compression
We construct the physical problem in a homogeneous ocean water environment where the standard density variance due to the water column's own weight is ignored, leading to ocean compressibility without static compression (only the dynamic compression is included). The mathematical problem is formulated in a two-dimensional coordinate system
$(x,z)$
whose
$x$
-direction is horizontal and
$z$
-direction vertically upwards. Linearised water wave theory is used to model the wave scattering induced by an initial pressure in a compressible ocean. The water is considered inviscid and the motion is irrotational, leading to the existence of a time-dependent velocity potential
$\varPhi (x,z,t)$
throughout the fluid domain such that
$\boldsymbol{\nabla }\varPhi (x,z,t)=(u(x,z,t),\,v(x,z,t)) =\boldsymbol{U}(x,z,t)$
, where
$u(x,z,t)$
and
$v(x,z,t)$
are the fluid velocities in the
$x$
and
$z$
directions, respectively. The equilibrium free surface and the rigid ocean floor are defined by
$z=0$
and
$z=-h$
, respectively.
Schematic diagram of the physical problem showing wave propagation in a compressible ocean (without including static compression) due to an underwater explosion (red star). Pressure waves generated by the explosion (red dashed lines) are reflected off the seabed and free surface (blue dashed lines) and induce propagating waves on the fluid's surface (solid red lines).
In the absence of extraneous force, the equations of motion can be written as
\begin{align} \boldsymbol{U}_t+(\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{U}&=-\frac {1}{\rho }\boldsymbol{\nabla }\!p=-\frac {1}{\rho } \frac {\partial p}{\partial \rho }\boldsymbol{\nabla }\!\rho , \\[-12pt]\nonumber \end{align}
\begin{align} \rho _t+\boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{U})&=0, \end{align}
where
$\rho$
and
$p$
represent the density and pressure, respectively. If
$s$
represents the relative change in the density from the undisturbed state (termed the ‘condensation’), one can write
$\rho =\rho _0(1+s)$
, where
$\rho _0$
represents the water density in the undisturbed state. Additionally, using the definition of the velocity potential and neglecting the higher-order terms, the above equations turn into
\begin{align} \boldsymbol{\nabla }\phi _t&=-\frac {1}{\rho _0} \frac {\text{d}p}{\text{d}\rho }\boldsymbol{\nabla }\!s, \\[-12pt]\nonumber \end{align}
\begin{align} s_t&=-\boldsymbol{\nabla} ^2\phi . \end{align}
Denoting
$c$
and
$K_0$
as sound speed and bulk modulus of the fluid, respectively, and following Lamb (Reference Lamb1957), we get
\begin{equation} \left .\frac {dp}{d\rho }\right |_{\rho =\rho _0} = \frac {K_0}{\rho _0}=c^2. \end{equation}
Eliminating
$s$
from (2.2) and making use of (2.3), we obtain the governing equation as
\begin{equation} \varPhi _{tt}=c^2\boldsymbol{\nabla} ^2\varPhi . \end{equation}
Hence, the IBVP, as given by Yamamoto (Reference Yamamoto1982) and Stiassnie (Reference Stiassnie2010), is
\begin{align} \varPhi _{tt}&=c^2\boldsymbol{\nabla} ^2\varPhi ,\quad &(x,z)\in \varOmega , \\[-12pt]\nonumber \end{align}
\begin{align} \varPhi _{tt}+g\varPhi _{z}&=0, \quad &z=0, \\[-12pt]\nonumber \end{align}
\begin{align} \varPhi _{z}&=0, &z=-h, \\[-12pt]\nonumber \end{align}
\begin{align} \varPhi _x&\to 0, & |x|\to \infty , \end{align}
subject to the initial conditions
\begin{align} \varPhi (x,z,0) &= \varPhi _0(x,z), \\[-12pt]\nonumber \end{align}
\begin{align} \varPhi _t(x,z,0)&=-\frac {P_0(x,z)}{\rho }, \end{align}
where
$c$
represents the speed of sound in sea water,
$g$
is acceleration due to gravity,
$\varOmega =\{(x,z): |x|\lt \infty ,-h\lt z\lt 0\}$
is the fluid domain and
$\varPhi _0$
and
$P_0$
are the initial potential and hydrodynamic pressure of the fluid, respectively.
We consider this problem in the context of self-adjoint operator theory. We rewrite (2.5a ) as
\begin{equation} \varPhi _{tt} - c^2{\nabla} ^2\varPhi =0, \end{equation}
and find an inner product with respect to which the operator
$-c^2\boldsymbol{\nabla} ^2$
is self-adjoint subject to the boundary conditions. To find this inner product, we use Green's second identity in the following form:
\begin{align} \iint _\varOmega \big (\boldsymbol{\nabla} ^2\varPhi \,\varPsi -\varPhi \,{\nabla} ^2\varPsi \big )\textrm {d}V&=\oint _{\partial \varOmega }\left (\varPhi _n\,\varPsi -\varPhi \,\varPsi _n\right )\textrm {d}S, \end{align}
where
$\varPhi$
and
$\varPsi$
are the two solutions of the IBVP (2.5a
)–(2.5d
). After rearranging the terms, we get
\begin{align} \iint _\varOmega \boldsymbol{\nabla} ^2\varPhi \,\varPsi \textrm {d}V-\oint _{\partial \varOmega }\varPhi _n\,\varPsi \textrm {d}S=\iint _\varOmega \varPhi \,\boldsymbol{\nabla} ^2\varPsi \textrm {d}V-\oint _{\partial \varOmega }\varPhi \,\varPsi _n\textrm {d}S. \end{align}
Taking
$\varOmega$
to be the region
$(-L,L)\times (-h,0)$
gives, in the limit as
$L\to \infty$
,
\begin{align} &\iint _\varOmega \boldsymbol{\nabla} ^2\varPhi \,\varPsi \textrm {d}V-\int _{-\infty }^\infty \varPhi _z\,\varPsi \,\textrm {d}x -\lim _{x\rightarrow \infty }\int _{-h}^0\varPhi _x\,\varPsi \,\textrm {d}z+\lim _{x\rightarrow -\infty }\int _{-h}^0\varPhi _x\,\varPsi \,\textrm {d}z \nonumber \\ &=\iint _\varOmega \boldsymbol{\nabla} ^2\varPsi \,\varPhi \textrm {d}V-\int _{-\infty }^\infty \varPsi _z\,\varPhi \,\textrm {d}x -\lim _{x\rightarrow \infty }\int _{-h}^0\varPsi _x\,\varPhi \,\textrm {d}z+\lim _{x\rightarrow -\infty }\int _{-h}^0\varPsi _x\,\varPhi \,\textrm {d}z, \end{align}
where we have used (2.5c ). Finally, use of (2.5d ) gives
\begin{align} \iint _\varOmega \boldsymbol{\nabla} ^2\varPhi \,\varPsi \textrm {d}V+\frac {c^2}{g}\int _{-\infty }^\infty \boldsymbol{\nabla} ^2\varPhi \,\varPsi \,\textrm {d}x=\iint _\varOmega \varPhi \,\boldsymbol{\nabla} ^2\varPsi \textrm {d}V+\frac {c^2}{g}\int _{-\infty }^\infty \varPhi \,\boldsymbol{\nabla} ^2\varPsi \,\textrm {d}x. \end{align}
Here, we used the condition on
$\varPhi$
and
$\varPsi$
and a combination of (2.5a
) and (2.5b
) that leads to
$\varPhi _z=-({c^2}/{g})\boldsymbol{\nabla} ^2\varPhi$
.
We define the Hilbert space to be the space equipped with the inner product
\begin{align} \langle \varPhi ,\varPsi \rangle _{\mathcal{E}}=\int _\varOmega \varPhi (x,z)\overline {\varPsi (x,z)}\,\textrm {d}x\textrm {d}z+\frac {c^2}{g}\int _{-\infty }^\infty \varPhi (x,0)\overline {\varPsi (x,0)}\,\textrm {d}x, \end{align}
for
$\varPhi$
,
$\varPsi \in \mathcal{E}$
. Self-adjoint operator theory requires an inner product space in which the operator has a symmetry, as we have just shown, and also that the domain of the operator and its adjoint are the same. The latter point is rather subtle, but in our case, this will follow when the function (or solution) space is closed. No formal proof of this is given. It is worth noting that the operator and the Hilbert space are closely related and that, when we close the function space, this closure occurs in the norm defined by our inner product. This function space is also the closure of the eigenfunctions of our operator, which we use to make our numerical calculations. In what follows, we use the inner product to calculate all projections, which is critical to the solution we develop.
3. Solution technique
Since the operator
$-c^2{\nabla} ^2$
is assumed to be self-adjoint, we know that its spectrum is real. In fact, it has a continuous spectrum occupying the entire positive real line. To find the generalised eigenfunctions, we solve
\begin{align} c^2\boldsymbol{\nabla} ^2\phi &= -\omega ^2\phi ,&(x,z)\in \varOmega , \\[-12pt]\nonumber \end{align}
\begin{align} -\omega ^2\phi +g\phi _z&=0,& z=0, \\[-12pt]\nonumber \end{align}
\begin{align} \phi _n&=0,& z=-h, \end{align}
subject to an incident wave condition. Solutions of (3.1a ) are nothing more than the frequency-domain solutions found by assuming time harmonicity, namely
\begin{equation} \varPhi (x,z,t) = {\textrm{Re}}\{\phi (x,z) \mathrm{e}^{-\textrm {i}\omega t}\}, \end{equation}
where
$\omega$
is the angular frequency.
Using separation of variables, we obtain explicit expressions for the solutions of (3.1a ) of the form
\begin{align} \phi ^{\pm }_n(x,z;k)=\mathrm{e}^{{\pm \textrm {i}}k x}f_n(z;k), \end{align}
for
$n\in \mathbb{N}$
and
$k\in \mathbb{R}$
, where the depth dependent functions satisfy the eigenvalue problem
\begin{equation} f_n''(z;k)=\mu _n^2f_n(z;k), \end{equation}
subject to the boundary conditions
\begin{align} gf_n'(0;k)&=c^2(k^2-\mu _n^2)f_n(0;k), \\[-12pt]\nonumber \end{align}
\begin{align} f_n'(-h;k)&=0, \end{align}
where the prime denotes differentiation with respect to
$z$
. Equation (3.4) has a discrete spectrum of eigenvalues
$\mu _n$
, indexed by
$n\in \mathbb{N}_0$
, satisfying the following dispersion relation:
\begin{align} \mu _n(k)\tanh {\mu _n(k) h}+\frac {c^2}{g^2}\mu _n(k)^2=\frac {c^2k^2}{g^2}, \end{align}
with
$\mu _0(k)$
being the positive real root and
$\mu _n(k)$
being of the form
${\textrm {i}}\tilde {\mu }_n(k)$
for
$n\geqslant 1$
, where
$\tilde {\mu }_n(k)$
are real, positive and in ascending order with respect to the index
$n$
. In addition,
$k_n$
are real positive for all progressive AGW modes and imaginary for evanescent modes. Note that, while
$\mu _n(k)$
depends on both
$n$
and
$k$
, we will omit the argument
$k$
when this is clear from context. The associated vertical eigenfunctions are written as
\begin{align} f_n(z;k)=&\frac {\displaystyle \cosh {\mu _n(z+h)}}{\cosh {(\mu _n h)}}. \end{align}
With reference to (3.1a
), the time harmonic angular frequency associated with the mode
$\phi _n(\boldsymbol{\cdot },\boldsymbol{\cdot };k)$
is
\begin{equation} \omega _n(k)^2= c^2(k^2-\mu _n^2). \end{equation}
Using the spectral theorem, the general solution of the time-domain problem (2.5) can be written as the following series of integrals over the continuous wavenumber:
\begin{equation} \varPhi (x,z,t)= \sum _{n=0}^{\infty }\int _{0}^\infty \left (C_n^-(t,k)\phi _n^{-}(x,z;k) \textrm {d}k + C_n^+(t,k) \phi _n^{+}(x,z;k)\right )\textrm {d}k. \end{equation}
Substitution of the above into (2.5a ) then gives
\begin{align} \sum _{n=0}^{\infty }\int _{0}^\infty -\omega _n(k)^2 \left (C_n^-(t,k)\phi _n^{-}(x,z;k) \textrm {d}k + C_n^+(t,k) \phi _n^{+}(x,z;k)\right )\textrm {d}k\nonumber \\ =\sum _{n=0}^{\infty }\int _{0}^\infty \left (C_n^{-\prime \prime }(t,k)\phi _n^{-}(x,z;k) \textrm {d}k + C_n^{+\prime \prime }(t,k) \phi _n^{+}(x,z;k)\right )\textrm {d}k, \end{align}
where the prime denotes differentiation with respect to
$t$
. Using the following orthogonality results (which we prove in Appendix A):
\begin{align} \langle \phi _n^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },\tilde {k})\rangle &= D_{\!p}(k)\delta (k-\tilde {k})\delta _{np}&k,\tilde {k}\gt 0,\quad n,p \in \mathbb{N}, \\[-12pt]\nonumber \end{align}
\begin{align} \langle \phi _n^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k),\phi _{\!p}^\mp (\boldsymbol{\cdot },\boldsymbol{\cdot },\tilde {k})\rangle &= 0&k,\tilde {k}\gt 0,\quad n,p \in \mathbb{N}, \end{align}
noting that (3.10a
) is the inner product between waves travelling in the same direction and (3.10b
) is the inner product between waves travelling in opposite directions and where
$D_{\!p}(k)$
are normalisation coefficients,
$\delta (\boldsymbol{\cdot })$
is the Dirac delta and
$\delta _{np}$
is the Kronecker delta, we obtain
\begin{equation} C_n^{\pm \prime \prime }(t,k)=-\omega _n(k)^2C_n^\pm (t,k). \end{equation}
Solutions to the above ordinary differential equation can be written as
\begin{equation} C_n^\pm (t,k)=A_n^\pm (k)\cos (\omega _n(k)t)+B_n^\pm (k)\frac {\sin (\omega _n(k)t)}{\omega _n(k)}, \end{equation}
whence
\begin{align} \varPhi (x,z,t) &= \sum _{n=0}^{\infty }\int _{0}^\infty \left (A_n^-(k)\phi _n^{-}(x,z;k)+A_n^+(k) \phi _n^{+}(x,z;k)\right ) \cos (\omega _n(k) t)\textrm {d}k \nonumber \\ &\quad + \sum _{n=0}^{\infty }\int _{0}^\infty \left (B_n^-(k)\phi _n^{-}(x,z;k)+B_n^+(k) \phi _n^{+}(x,z;k)\right ) \frac {\sin (\omega _n(k) t)}{\omega _n(k)}\textrm {d}k. \end{align}
To obtain the particular solution to the IBVP (2.5), the coefficient functions
$A_n^\pm$
and
$B_n^\pm$
must be determined from the initial potential and pressure of the fluid,
$\varPhi _0$
and
$P_0$
, respectively. Applying (3.10)–(3.13) for
$t=0$
yields
\begin{align} \langle \varPhi (\boldsymbol{\cdot },\boldsymbol{\cdot },0),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot }, \tilde {k}) \rangle _{\mathcal{E}} &=\sum _{n=0}^{\infty } \int _{0}^\infty A_n^-(k)\langle \phi _n^-(\boldsymbol{\cdot },\boldsymbol{\cdot },k), \phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },\tilde {k})\rangle _{\mathcal{E}} \nonumber \\&\quad +A_n^+(k)\langle \phi _n^-(\boldsymbol{\cdot },\boldsymbol{\cdot },k),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },\tilde {k})\rangle _{\mathcal{E}} \,\textrm {d}k, \nonumber \\ &=D_{\!p}(\tilde {k})A_{\!p}^\pm (\tilde {k}). \end{align}
Similarly
\begin{equation} \langle \varPhi _t(\boldsymbol{\cdot },\boldsymbol{\cdot },0),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },\tilde {k})\rangle _{\mathcal{E}} =D_{\!p}(\tilde {k})B_{\!p}^\pm (\tilde {k}). \end{equation}
Now we apply the initial conditions (2.5e
) and (2.5f
) and take inner products with respect to
$\phi _{\!p}^\pm$
to get
\begin{align} \langle \varPhi _0(x,z),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k)\rangle _{\mathcal{E}} &=\langle \varPhi (\boldsymbol{\cdot },\boldsymbol{\cdot },0),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k) \rangle _{\mathcal{E}}=D_{\!p}A_{\!p}^\pm (k), \\[-12pt]\nonumber \end{align}
\begin{align} -\frac {\langle P_0(x,z),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k)\rangle _{\mathcal{E}} }{\rho }&=\langle \varPhi _t(\boldsymbol{\cdot },\boldsymbol{\cdot },0),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k)\rangle _{\mathcal{E}}=D_{\!p}B_{\!p}^\pm (k). \end{align}
Then the coefficient functions
$A_{\!p}^\pm (k)$
and
$B_{\!p}^\pm (k)$
are given by
\begin{equation} A_{\!p}^\pm (k)=\frac {\langle \varPhi _0(x,z),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k)\rangle _{\mathcal{E}} }{D_{\!p}} \quad \mbox{and}\quad B_{\!p}^\pm (k)=-\frac {\langle P_0(x,z),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k)\rangle _{\mathcal{E}} }{\displaystyle \rho D_{\!p}}. \end{equation}
Hence, we obtain the potential function as
\begin{align} \varPhi (x,z,t)=\sum _{n=0}^{\infty }\int _{0}^\infty \left (\frac {\langle \varPhi _0(x,z),\phi _n^-\rangle _{\mathcal{E}}}{D_{\!p}} \phi _n^-+\frac {\langle \varPhi _0(x,z),\phi _n^+\rangle _{\mathcal{E}}}{D_{\!p}} \phi _n^+\right )\cos (\omega _n t)\,\textrm {d}k \nonumber \\ -\frac {1}{\rho }\sum _{n=0}^{\infty }\int _{0}^\infty \left (\frac {\langle P_0(x,z),\phi _n^-\rangle _{\mathcal{E}}}{D_{\!p}} \phi _n^-+\frac {\langle P_0(x,z),\phi _n^+\rangle _{\mathcal{E}}}{D_{\!p}} \phi _n^+\right )\frac {\sin (\omega _n t)}{\omega _n}\,\textrm {d}k. \end{align}
We note at this point some special features of this solution for
$\varPhi$
. The solution can be calculated at any time from the initial conditions without reference to the solution at any other time. Moreover, the difficulty of solving for any time is equivalent as the method does not rely on temporal iteration. In particular, we compute the solution at
$t=0$
, which provides a validation of the method through comparison with the initial condition. We note that, provided that the
$\phi _n^{\pm }$
are correct, then
$\varPhi$
must satisfy the governing equations, and if the potential matches at
$t=0$
, we have calculated the solution correctly.
4. Time-domain simulation
We will now compute time-domain solutions of the internal pressure
\begin{equation} P(x,z,t)=-\rho \varPhi _t(x,z,t), \end{equation}
and free surface
\begin{equation} \eta (x,t)=-\frac {1}{g}\varPhi _t(x,0,t)=\frac {1}{\rho g}P(x,0,t). \end{equation}
Noting that when computing the pressure over the full
$z$
-domain, we recover the surface as a byproduct.
4.1. Numerical evaluation
When evaluating (4.1), care must be taken due to the potentially highly oscillatory nature of the integrand (for large
$x$
or
$t$
). We reformulate (3.17) into the form
\begin{align} \varPhi (x,z,t) = \sum _{n=0}^\infty I_n^{(1)}-\frac {1}{\rho }\sum _{n=0}^\infty I_n^{(2)}, \end{align}
where
\begin{align} I_n^{(1)}&=\sum _{\pm }\sum _{j=\{0,1\}}(-1)^j\int _{0}^\infty g_{n\pm }^{(1)}{e}^{\textrm {i}(\pm kx+(-1)^j\omega _nt)}\,\textrm {d}k , \\[-12pt]\nonumber \end{align}
\begin{align} I_n^{(2)}&=\sum _{\pm }\sum _{j=\{0,1\}}\int _{0}^\infty g_{n\pm }^{(2)}{e}^{\textrm {i}(\pm kx+(-1)^j\omega _nt)}\,\textrm {d}k , \end{align}
and
\begin{align} g_{n\pm }^{(1)}=\frac {\omega _n\langle \varPhi _0(x,z),\phi _n^\pm \rangle _{\mathcal{E}}f_n(z)}{2D_{\!p}}, \quad \quad g_{n\pm }^{(2)}=\frac {\langle P_0(x,z),\phi _n^\pm \rangle _{\mathcal{E}}f_n(z)}{2D_{\!p}}. \end{align}
The infinite sum is truncated at
$N$
modes while the semi-infinite integral is truncated at
$k=k_{\textit{max}}$
and discretised uniformly using spacing
$\Delta k$
. The choice of
$N$
,
$k_{\textit{max}}$
and
$\Delta k$
depends on the initial conditions and will be discussed in § 4.2. In this paper, we take
$N=100$
as that was sufficient to reproduce the initial condition. Similarly, the semi-infinite integral is truncated at
$k=0.2$
and discretised uniformly using
$1001$
points.
This formulation allows us to implement a Filon-type quadrature method (Filon Reference Filon1930; Dahlquist & Björck Reference Dahlquist and Björck2008) detailed in Appendix B. Further efficiency can be gained by recognising that
\begin{align} g_{+}{e}^{\textrm {i}(kx+\omega _nt)}=\overline {g_{-}{e}^{-\textrm {i}(kx+\omega _nt)}}, \end{align}
halving the number of terms in the summations that must be computed.
4.2. Underwater pressure disturbance
As an example, we will consider an underwater pressure disturbance as the initial condition (generated by an underwater explosion, for example). The fluid is assumed initially at rest,
$\varPhi _0(x,z)=0$
, with the initial pressure distribution
\begin{equation} P_0(x,z) = \hat {P}e^{\displaystyle -\pi ^2\frac {(x-x_c)^2+(z-z_c)^2}{\sigma ^2}}, \end{equation}
where
$\hat {P}$
is the maximum pressure,
$(x_c,z_c)$
is the centre of the pressure distribution and
$\sigma$
is the characteristic width of the pressure. For the following results, we choose
$\hat {P}=10^6$
Pa,
$(x_c,z_c)=(0,-2000)$
m,
$g=9.81$
ms
$^{-2}$
,
$c=1450$
ms−1,
$h=4000$
m and
$\sigma =200$
m.
The numerical parameters
$N$
,
$k_{\textit{max}}$
and
$\Delta k$
are chosen by ensuring the initial condition can be reproduced with sufficient accuracy. We observe the maximum difference between the numerical solution and true initial condition over different parameter choices (see Appendix C for details). Based on this parameter sweep, we choose
$N=100$
,
$k_{\textit{max}}=0.2$
and
$\Delta k=0.0002$
.
Compression waves caused by an initial Gaussian pressure distribution centred at
$(0,-2000)$
.
Figure 2 shows the solutions at
$t=1$
,
$2$
and
$6.75$
s. We can see that the initial pressure spreads radially away from the initial location at the speed of sound. The expanding wave reflects off both the bottom and the free surface, the phase of the pressure wave reversing as it reflects off the surface (figure 2
b). This process continues as time evolves, resulting in multiple pressure waves propagating away from the pressure origin. The effect this has on the free surface is the creation of wave crests above the initial pressure location that propagate away. The free-surface elevation between the crests will also be greater than at equilibrium. Full animation is given as Movie 1.
5. Inclusion of static compression
In this section, we include (mild) density variance of the ocean water, the effect which is known in the literature as static compression, into the mathematical formulation, and develop the solution using a similar approach to that used in the previous section. The rigorous derivation of the governing equations can be found in the works of Longuet-Higgins (Reference Longuet-Higgins1950), Miyoshi (Reference Miyoshi1954) and Sells (Reference Sells1965). Note that, in the subsequent mathematical treatment, the density
$\rho$
represents its value at the surface, which is treated as the baseline density. The hydrostatic density is
$\rho e^{-\gamma z}$
(Abdolali & Kirby Reference Abdolali and Kirby2017), where the effect of the variance in the density away from the surface is reflected through the term
$\gamma$
, which depends on the speed of sound in the water and is defined by
$\gamma = g/c^2$
. This inclusion inherently modifies the underlying governing equation, which we briefly show now.
Assuming
$\rho _0(z)$
and
$p_0(z)$
being the undisturbed water density and pressure, respectively, the fundamental equations of motion up to the first order yield (Bondi Reference Bondi1947)
\begin{align} &\phi _t=-\int \frac {dp}{\rho }+gz, \\[-12pt]\nonumber \end{align}
\begin{align} &\rho _t+\boldsymbol{\nabla }\boldsymbol{\cdot }(\rho _0\boldsymbol{\nabla }\phi )=0. \end{align}
According to Longuet-Higgins (Reference Longuet-Higgins1950), pressure and density are related by the relationship
\begin{equation} \frac {\partial p}{\partial \rho }= c^2, \end{equation}
which modifies (5.1a ) to
\begin{equation} \phi _t=-c^2\log \rho +gz. \end{equation}
For the undisturbed state, one can obtain from (5.3)
\begin{equation} \rho _0= \hat {\rho }_0e^{gz/c^2} =\hat {\rho }_0e^{\gamma z}, \end{equation}
where
$\hat {\rho }_0$
represents undisturbed density at
$z=0$
. Using
$\rho =\rho _0(1+s)$
, (5.3) and (5.1b
) become
\begin{align} &\frac {\partial \phi }{\partial t}=-c^2s, \\[-12pt]\nonumber \end{align}
\begin{align} & \frac {\partial s}{\partial t}+\boldsymbol{\nabla} ^2\phi + \gamma \phi _z=0. \end{align}
Eliminating
$s$
yields the governing equation
\begin{equation} \boldsymbol{\nabla} ^2\varPhi = \frac {1}{c^2}\varPhi _{tt}+\gamma \varPhi _{z}. \end{equation}
The effect of static compression is small for realistic oceans (Abdolali et al Reference Abdolali, Kadri and Kirby2019), as we shall see, and it significantly complicates the formulation. We believe that including it is worthwhile, first because it allows accurate calculation of its effect and second because it highlights how powerful the special inner product-based method is and how it can be extended to highly non-trivial cases. We note that the complexity arises in the theory, but that the computational difficulty is not greatly increased.
5.1. Mathematical solution with static compression
The modified IBVP is given by
\begin{align} \boldsymbol{\nabla} ^2\varPhi &= \frac {1}{c^2}\varPhi _{tt}+\gamma \varPhi _{z},\quad (x,z)\in \varOmega , \\[-12pt]\nonumber \end{align}
\begin{align} \varPhi _{tt}+g\varPhi _{z}&=0, \qquad \qquad \qquad \qquad z=0, \\[-12pt]\nonumber \end{align}
\begin{align} \varPhi _{z}&=0, \qquad \qquad \qquad \qquad z=-h, \\[-12pt]\nonumber \end{align}
\begin{align} \varPhi _{x}&\to 0, \qquad \qquad \qquad \qquad |x|\to \infty , \end{align}
subject to the initial conditions
\begin{align} \varPhi (x,z,0) &= 0, \\[-12pt]\nonumber \end{align}
\begin{align} \varPhi _t(x,z,0)&=-\frac {P_0(x,z)}{\rho e^{-\gamma z}}. \end{align}
In order to identify an appropriate inner product, we apply Green's second identity to the functions
$\hat {\varPhi }=e^{-\gamma z/2}\varPhi$
and
$\hat {\varPsi }=e^{-\gamma z/2}\varPsi$
, where
$\varPhi$
and
$\varPsi$
satisfy (5.7), in the form
\begin{align} \iint _\varOmega \big (\boldsymbol{\nabla} ^2\hat {\varPhi }\,\hat {\varPsi }-\hat {\varPhi }\,\boldsymbol{\nabla} ^2 \hat {\varPsi }\big )\,\textrm {d}V=\oint _{\partial \varOmega }\big (\hat {\varPhi }_n\, \hat {\varPsi }-\hat {\varPhi }\,\hat {\varPsi }_n\big )\,\textrm {d}S. \end{align}
The choices of
$\hat {\varPhi }$
and
$\hat {\varPsi }$
ensure elimination of the first-order derivative term in (5.7a
). Now taking
$\varOmega$
to be the region
$(-L,L)\times (-h,0)$
gives, in the limit as
$L\to \infty$
\begin{align} \iint _\varOmega e^{-\gamma z}\big [\big (\boldsymbol{\nabla} ^2\varPhi -\gamma \varPhi _z\big )\varPsi -\varPhi \big (\boldsymbol{\nabla} ^2\varPsi -\gamma \varPsi _z\big )\big ]\,\textrm {d}V = \int _{-\infty }^\infty \left .\left (\varPhi _z\,\varPsi -\varPhi \,\varPsi _z\right )\right |_{z=0}\textrm {d}x, \end{align}
where we have used (5.7c ) and (5.7d ). Subsequently
\begin{align} \iint _\varOmega e^{-\gamma z}\big (\boldsymbol{\nabla} ^2\varPhi -\gamma \varPhi _z\big )\varPsi \textrm {d}V-\int _{-\infty }^\infty \left .\varPhi _z\,\varPsi \right |_{z=0}\textrm {d}x\;&\nonumber \\ =\iint _\varOmega e^{-\gamma z}\varPhi \big (\boldsymbol{\nabla} ^2\varPsi -\gamma \varPsi _z\big )\,\textrm {d}V-\int _{-\infty }^\infty \left .\varPhi \,\varPsi _z\right |_{z=0}\textrm {d}x.& \end{align}
Note that the combination of (5.7a
) and (5.7b
) yields
$\varPhi _z=-({c^2}/{g})(\boldsymbol{\nabla} ^2\varPhi -\gamma \varPhi _z)$
on
$z=0$
, which further modifies the above expression into
\begin{align} \iint _\varOmega e^{-\gamma z}\mathcal{A}\varPhi \varPsi \textrm {d}V+\frac {c^2}{g}\int _{-\infty }^\infty \left .\mathcal{A}\varPhi \,\varPsi \right |_{z=0}\textrm {d}x&=\iint _\varOmega e^{-\gamma z}\varPhi \mathcal{A}\varPsi \textrm {d}V+\frac {c^2}{g}\int _{-\infty }^\infty \left .\varPhi \,\mathcal{A}\varPsi \right |_{z=0}\textrm {d}x, \end{align}
where
$\mathcal{A}=\boldsymbol{\nabla} ^2-\gamma \partial _z$
. Consequently, this operator is symmetric provided we choose the inner product to be
\begin{equation} \langle \varPhi ,\varPsi \rangle _{\mathcal{E}} =\int _\varOmega e^{-\gamma z}\varPhi (x,z)\overline {\varPsi (x,z)}\,\textrm {d}x\textrm {d}y+\frac {c^2}{g}\int _{-\infty }^\infty \varPhi (x,0)\overline {\varPsi (x,0)}\,\textrm {d}x. \end{equation}
In a way similar to what is mentioned in § 3, we expand the time-domain solution in the generalised eigenfunctions of the underlying differential operator. These eigenfunctions satisfy
\begin{align} \boldsymbol{\nabla} ^2\phi &= -\frac {\omega ^2}{c^2}\phi +\gamma \phi _{z},\quad &(x,z)\in \varOmega , \\[-12pt]\nonumber \end{align}
\begin{align} -\omega ^2\phi +g\phi _{z}&=0, \quad &z=0, \\[-12pt]\nonumber \end{align}
\begin{align} \phi _{z}&=0, &z=-h, \end{align}
and are of the form
\begin{equation} \phi _n^\pm (x,z;k)={e}^{\textrm {i}kx}f_n(z;k). \end{equation}
In the case of static compression, the depth-dependent functions
$f_n(z,k)$
are more complicated. They satisfy the eigenvalue problem
\begin{align} f_n''(z;k)-\gamma f_n'(z;k) &= \left (\mu _n^2-\frac {\gamma ^2}{4}\right )f_n(z;k), \\[-12pt]\nonumber \end{align}
\begin{align} f_n'(-h;k)&=0, \\[-12pt]\nonumber \end{align}
\begin{align} f'(0;k)&=\frac {c^2}{g}\left (k^2+\frac {\gamma ^2}{4}-\mu _n^2\right )f(0;k), \end{align}
where the prime denotes differentiation with respect to
$z$
. For each
$k$
, the discrete spectrum of eigenvalues
$\mu _n$
of (5.16) satisfy the following dispersion relation:
\begin{align} \frac {c^2}{g}\left (k^2+\frac {\gamma ^2}{4}-\mu _n^2\right )\left (\displaystyle 1-\left (\frac { \gamma }{2\mu _n} \right )\tanh {\mu _n h}\right )=\mu _n\! \displaystyle \left [1- \left ( \frac {\gamma }{2\mu _n} \right )^2 \right ]\tanh {\mu _n h}. \end{align}
The solutions to (5.17) are taken so that
$\mu _0(k)\in \mathbb{R}$
and
$\mu _n(k)$
are purely imaginary with positive and increasing imaginary part. The angular frequency associated with the mode
$\phi _n^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot };k)$
is
\begin{equation} \omega _n(k)=c\sqrt {k^2+\frac {\gamma ^2}{4}-\mu _n^2.} \end{equation}
Without loss of generality, we set
$f(0)=1$
, yielding
\begin{align} f_n(z;k)=&\frac {\displaystyle e^{(\gamma z/2)}\left ( \frac {\gamma }{2\mu _n}\sinh {\mu _n(z+h)}-\cosh {\mu _n(z+h)} \right )}{\dfrac {\gamma }{2\mu _n}\sinh {( \mu _n h)}-\cosh {(\mu _n h)}}. \end{align}
As in § 3, the spectral theorem gives a time-domain solution of the form
\begin{align} \varPhi (x,z,t)&= \sum _{n=0}^{\infty }\int _{0}^\infty \left (A_n^-(k)\phi _n^{-}(x,z;k)+A_n^+(k) \phi _n^{+}(x,z;k)\right ) \cos (\omega _n(k) t)\textrm {d}k \nonumber \\ &\quad + \sum _{n=0}^{\infty }\int _{0}^\infty \left (B_n^-(k)\phi _n^{-}(x,z;k)+B_n^+(k) \phi _n^{+}(x,z;k)\right ) \frac {\sin (\omega _n(k) t)}{\omega _n(k)}\textrm {d}k, \end{align}
where the coefficient functions
$A_n^\pm$
and
$B_n^\pm$
are to be determined from the initial conditions. With orthonormality conditions again being of the form (3.10), albeit with different normalisation functions
$D_{\!p}(k)$
, the coefficients are recovered as
\begin{align} A_{\!p}^\pm (k)=\frac {\langle \varPhi _0(x,z),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k)\rangle _{\mathcal{E}} }{D_{\!p}} \quad \mbox{and}\quad B_{\!p}^\pm (k)=-\frac {\langle P_0(x,z),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k)\rangle _{\mathcal{E}} }{\displaystyle D_{\!p}}, \end{align}
which is precisely in the form of (3.16). To finalise the mathematical solution in the case of static compression, it remains to compute the normalisation functions
$D_{\!p}(k)$
. This is done in Appendix D. We note that (5.20) and (5.21) reduce to (3.17) and (3.16), respectively, by setting
$\gamma =0$
.
5.2. Time-domain simulations with static compression
Compression waves caused by an initial Gaussian pressure distribution centred at
$(0,-2000)$
, including the effects of static compression. In each panel, the line graph represents the free-surface displacement, and the colour plot represents the hydrodynamic pressure distribution below the surface.
The maximum percentage pressure difference (5.22) for different depths,
$h$
, where the pressure is centred either half-way between the bottom and the surface
$(x_c,z_c)=(0,-h/2)$
, or at
$(0,-250)$
.
A plot of the pressure difference between figures 2(a) and 3(a) as a percentage of that maximum magnitude of pressure at that time (
$t=1$
), where a positive difference indicates that the pressure solution including static compression is greater than the pressure solution without static compression. The bottom panel provides a cross-section of the pressure difference at
$x=0$
.
Firstly, we consider the underwater explosion example detailed in § 4.2, to compare with the results in the absence of static compression. We can see in figure 3 (the full animation is in Movie 2) that the results are visually indistinguishable from the prior results without static compression (figure 2). To better observe the effects of static compression, we plot the percentage difference
\begin{equation} \frac {P_{\textit{static}}(x,z,t)-P(x,z,t)}{\max \limits _{x,z}(P_{\textit{static}}(x,z,t),P(x,z,t))}, \end{equation}
for
$t=1$
(figure 4), where
$P_{\textit{static}}$
is the solution including static compression. We observe that the difference wave has a negative leading front in the upper half-wave (i.e. upwards travelling waves) and a positive leading front in the lower half (i.e. downward travelling waves). We can see that the inclusion of static compression, on average, causes the magnitude of pressure above (below) the initial explosion location to decrease (increase) compared with the case with no static compression. It aligns with the fact that water is less (more) dense in the upper (lower) water region measured from the point of explosion when static compression is considered. Further, we calculated the maximum of the percentage difference (5.22) over space (
$|x|\lt 10^5$
,
$-h\lt z\lt 0$
) and time (
$0\lt t\lt 20$
) for different depths,
$h$
, where the pressure is centred half-way between the bottom and the surface
$(x_c,z_c)=(0,-h/2)$
and for a pressure centre fixed at
$(0,-250)$
. From the results presented in table 1, we can see that the pressure difference increases as the depth decreases, and in the case of the pressure at
$z_c=-250$
m, there is a minimum in pressure difference due to static compression at roughly
$h=1500$
m before the difference increases as the depth increases. However, the differences still remain roughly within
$\pm 1\,\%$
for all tested depths. Finally, we note that the pressure wave speed has not changed with the addition of static compression. The Movie 3 provides the full animation for figure 4.
Compression waves caused by an initial Gaussian pressure distribution centred at
$(0,0)$
, including the effects of static compression.
We will now consider two different initial pressure distributions that have the same surface displacement but drastically different subsurface profiles. The first initial pressure distribution is the Gaussian distribution (4.8) centred at the free-surface
$(x_c,z_c)=(0,0)$
, and the second is a one-dimensional Gaussian with the same
$x$
behaviour as (4.8)
\begin{equation} P_1(x,z) = \hat {P}\exp\! \left (\!{\displaystyle -\pi ^2\frac {(x-x_c)^2}{\sigma ^2}}\right )\!. \end{equation}
Comparing the results for the localised Gaussian (figure 5) with those for the one-dimensional Gaussian (figure 6), we can see that the initial pressures excite propagating gravity waves on the surface of the same shape, which move much slower than the pressure waves. These gravity waves are orders of magnitude greater than the surface displacement due to compression waves. This means that, in the case that gravity waves and pressure waves are created by the same source, it would be hard to detect the effects of compression and further determine the initial pressure profile from surface observations. The pressure due to the gravity waves decays away from the surface, and therefore, underwater pressure measurements are better suited to observe the internal compression waves.
Focusing now on the internal pressure waves for our two initial pressures, as can be expected, the radically different initial pressure distributions lead to different internal pressure behaviours. The pressure wave generated by the localised Gaussian (figure 5) propagates radially, reflecting off the bottom and surface, creating a similar surface knocking phenomenon to the submerged Gaussian. The one-dimensional Gaussian (figure 6), on the other hand, predominately propagates in the
$x$
-direction with a rapidly decaying trail that reflects off the surface, which does not produce an observable knocking behaviour. The respective animations are provided in Movie 4 and Movie 5.
Compression waves caused by an initial one-dimensional Gaussian pressure distribution centred at
$x_c=0$
, including the effects of static compression.
6. Conclusion
Under the assumptions of linearised water wave theory, this paper has considered the IBVP governing the evolution of a pressure disturbance in a two-dimensional ocean with a free surface, including the effect of compression. The problem was solved mathematically using self-adjoint operator theory by introducing an appropriate inner product. We presented two versions of the theory; in the first, we ignored the effect of static compression and included only the dynamic compression, while the second included both dynamic and static compression. Time-domain simulations showed the propagation of the pressure wave away from the initial source, reflecting off the free surface and the rigid ocean floor, and the onset of surface gravity waves propagating away from the source. In particular, a near reversal of phase is observed upon reflection from the surface, while the opposite pattern occurs in the reflection process from the rigid ocean floor. This indicates the fact that the free surface is approximately a homogeneous Dirichlet boundary. The corresponding time-domain simulations for the model incorporating static compression show that the effect of static compression is minor. While the present work deals with the canonical problem of uniform water depth, the theory considered here can be extended to tackle more realistic variable bathymetry using step approximation, and it is equally applicable in modelling three-dimensional waves, albeit with increased computational requirements.
Supplementary material and movies
Animations for figures 2, 3, 4, 5 and 6 are uploaded as supplementary material. Supplementary movies are available at https://doi.org/10.1017/jfm.2026.11386.
Funding
B.W. and M.H.M. are supported by the Australian Research Council's Discovery Projects funding scheme (DP240102104).
Data availability
Data sharing does not apply to this article, as no new data were created or analysed in this study. The computer code used to generate the solutions is available from the corresponding author upon reasonable request.
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Orthogonality of the eigenfunctions
Here, we determine the normalisation coefficients
$D_{\!p}(k)$
while simultaneously confirming the orthogonality of the eigenfunctions. Firstly, we calculate the
$L^2(-h,0)$
inner product of the depth-dependent functions through two applications of integration by parts, yielding
\begin{align} \int _{-h}^0f_n(z;k)f_{\!p}(z;\tilde {k})\,\textrm {d}z&=\int _{-h}^0 \frac {\displaystyle \cosh {\mu _n(z+h)}}{\cosh {\mu _n h}}\frac {\displaystyle \cosh {\tilde {\mu }_{\!p}(z+h)}}{\cosh {\tilde {\mu }_{\!p} h}}\,\textrm {d}z\nonumber \\ & =\frac {1}{\mu _n}\tanh {\mu _n h}-\frac {\tilde {\mu }_{\!p}}{\mu _n^2}\tanh {\tilde {\mu }_{\!p} h}+\frac {\tilde {\mu }_{\!p}^2}{\mu _n^2}\int _{-h}^0f_n(z;k)f_{\!p}(z;\tilde {k})\,\textrm {d}z, \end{align}
where for shorthand we write
$\mu _n=\mu _n(k)$
and
$\tilde {\mu _{\!p}}=\mu _{\!p}(\tilde {k})$
. Rearranging using (3.5) and (3.7) to remove instances of
$\mathrm{tanh}(\boldsymbol{\cdot })$
and
$k$
gives
\begin{align} \int _{-h}^0f_n(z;k)f_{\!p}(z;\tilde {k})\,\textrm {d}z=\frac {1}{g}\frac {\omega _n^2-\tilde {\omega }_{\!p}^2}{\mu _n^2-\tilde {\mu }_{\!p}^2}, \end{align}
with
$\omega _n=\omega _n(k)$
and
$\tilde {\omega _{\!p}}=\omega _{\!p}(\tilde {k})$
. Now we compute the inner product
$\langle \phi _n^\pm ,\phi _{\!p}^\pm \rangle _{\mathcal{E}}$
as
\begin{align} \langle \phi _n^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },\tilde {k}) \rangle _{\mathcal{E}} &=\int _{-\infty }^\infty \int _{-h}^0 \varphi _n^\pm (x,z;k)\overline {\varphi _{\!p}^\pm (x,z;\tilde {k})}\,\textrm {d}z\,\textrm {d}x\nonumber \\ &\quad \quad +\frac {c^2}{g}\int _{-\infty }^\infty \varphi _n^\pm (x,0;k) \overline {\varphi _{\!p}^\pm (x,0;\tilde {k})}\,\textrm {d}x,\nonumber \\ &=\int _{-\infty }^\infty \left (\frac {\omega _n^2-\tilde {\omega }_{\!p}^2}{\mu _n^2-\tilde {\mu }_{\!p}^2}+\frac {c^2}{g}\right )e^{\mp \textrm {i}(k-\tilde {k}) x}\,\textrm {d}x,\nonumber \\ &=\frac {c^2}{g}\int _{-\infty }^\infty \left (\frac {k^2-\tilde {k}^2}{\mu _n^2-\tilde {\mu }_{\!p}^2}\right )e^{\mp \textrm {i}(k-\tilde {k}) x}\,\textrm {d}x,\nonumber \\ &=2\pi \frac {c^2}{g}\left (\frac {k^2-\tilde {k}^2}{\mu _n^2-\tilde {\mu }_{\!p}^2}\right )\delta (k-\tilde {k}), \end{align}
where
$\delta (\kappa )$
is the Dirac delta. Finally, in the case of
$n=p$
, we can recover the normalisation constant by taking limits
\begin{align} \lim _{k\rightarrow \tilde {k}}\langle \phi _n^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },\tilde {k})\rangle _{\mathcal{E}}= 2\pi \frac {c^2}{g}\left (\frac {\tilde {k}}{\mu _n\partial _k\tilde {\mu }_n}\right )\delta (k-\tilde {k})\delta _{np}, \end{align}
where
\begin{align} \partial _k \mu _n=\frac {{2}\frac{c^2}{g}k}{\mu _n h\textrm {sech}^2{\mu _n h}+\tanh {\mu _n h}+2\dfrac {c^2}{g}\mu _n}. \end{align}
For simplicity, we write
\begin{align} \lim _{k\rightarrow \tilde {k}}\langle \phi _n^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },\tilde {k})\rangle _{\mathcal{E}}= D_{\!p}(k)\delta (k-\tilde {k})\delta _{np}, \end{align}
where
\begin{align} D_{\!p}(k)=2\pi \left (\frac {h}{2}\textrm {sech}^2{\mu _{\!p} h}+\frac {\tanh {\mu _{\!p} h}}{2\mu _{\!p}}+\frac {c^2}{g}\right ). \end{align}
Similarly,
\begin{align} \langle \phi _n^\pm ,\phi _{\!p}^\mp \rangle _{\mathcal{E}} =2\pi \frac {c^2}{g}\left (\frac {k^2-\tilde {k}^2}{\mu _n^2-\tilde {\mu }_{\!p}^2}\right )\delta (k+\tilde {k})=0, \end{align}
since
$k$
and
$\tilde {k}$
are both positive real numbers.
Appendix B. Filon-type quadrature
Following the method of Filon (Reference Filon1930), we compute an integral of the form
\begin{align} I=\int _{a}^{b}f(x){e}^{\textrm {i}\psi (x)}\,\text{d}x, \end{align}
where
$f(x)$
and
$\psi '(x)$
are slowly varying functions, by splitting the domain of integration at the points
$x_0=a,x_1,\ldots ,x_N=b$
and approximating the functions
$f(x)$
and
$\psi (x)$
by piecewise linear functions
\begin{align} I&=\sum _{i=0}^{N-1}\int _{x_i}^{x_{i+1}}\left (f_i+\frac {f_{i+1}-f_i}{\Delta x_i}(x-x_i)\right ){e}^{\textrm {i}\left (\psi _i+\frac {\psi _{i+1}-\psi _i}{\Delta x_i}(x-x_i)\right )}\text{d}x,\nonumber \\ &=\sum _{i=0}^{N-1}\Delta x_i\int _{0}^{1}\left (f_i+(f_{i+1}-f_i)s\right ){e}^{\textrm {i}\left (\psi _i+(\psi _{i+1}-\psi _i)s\right )}ds,\quad \mbox{where } \quad s=\frac {x-x_i}{\Delta x_i},\nonumber \\ &=\sum _{i=0}^{N-1}\Delta x_i\left [-\frac {\textrm {i}}{\psi _{i+1}-\psi _i}\big (f_{i+1}{e}^{\textrm {i}\psi _{i+1}} -f_i{e}^{\textrm {i}\psi _i}\big ) +\frac {f_{i+1}-f_i}{(\psi _{i+1}-\psi _i)^2}\big ({e}^{\textrm {i}\psi _{i+1}}-{e}^{\textrm {i}\psi _i}\big )\right ]\!,\nonumber \\ &=\sum _{i=0}^{N-1}\Delta x_i\left [\left (\frac {\textrm {i}{e}^{\textrm {i}\psi _i}}{\psi _{i+1}-\psi _i}-\frac {{e}^{\textrm {i}\psi _{i+1}}-{e}^{\textrm {i}\psi _i}}{(\psi _{i+1}-\psi _i)^2}\right )f_i +\left (-\frac {\textrm {i}{e}^{\textrm {i}\psi _{i+1}}}{\psi _{i+1}-\psi _i}+\frac {{e}^{\textrm {i}\psi _{i+1}}-{e}^{\textrm {i}\psi _i}}{(\psi _{i+1}-\psi _i)^2}\right )f_{i+1}\right ]\!,\nonumber \\ &=\Delta x_0\left (\frac {\textrm {i}{e}^{\textrm {i}\psi _0}}{\psi _{1}-\psi _0}-\frac {{e}^{\textrm {i}\psi _{1}}-{e}^{\textrm {i}\psi _0}}{(\psi _{1}-\psi _0)^2}\right )f_0+\sum _{i=1}^{N-1}\left [\Delta x_{i-1}\left (-\frac {\textrm {i}{e}^{\textrm {i}\psi _{i}}}{\psi _{i}-\psi _{i-1}}+\frac {{e}^{\textrm {i}\psi _{i}}-{e}^{\textrm {i}\psi _{i-1}}}{(\psi _{i}-\psi _{i-1})^2}\right )\right .\nonumber \\ &\quad \left .+\Delta x_i\left (\frac {\textrm {i}{e}^{\textrm {i}\psi _i}}{\psi _{i+1}-\psi _i}-\frac {{e}^{\textrm {i}\psi _{i+1}}-{e}^{\textrm {i}\psi _i}}{(\psi _{i+1}-\psi _i)^2}\right )\right ]f_i+\Delta x_{N-1}\left (-\frac {\textrm {i}{e}^{\textrm {i}\psi _{N}}}{\psi _{N}-\psi _{N-1}}+\frac {{e}^{\textrm {i}\psi _{N}}-{e}^{\textrm {i}\psi _{N-1}}}{(\psi _{N}-\psi _{N-1})^2}\right )f_{N}. \end{align}
Initial condition reconstruction error for different
$N$
, and
$k_{\textit{max}}=0.2$
and
$\Delta k=0.0002$
.
The result is a weighting scheme for quadrature integration that can account for a rapidly changing phase function.
Appendix C. Errors for numerical parameter choice
This section compares the effects that the numerical scheme parameters
$N$
,
$k_{\textit{max}}$
and
$\Delta k$
have on the error in reconstructing the initial condition for the underwater Gaussian pressure disturbance (4.8). The physical parameters used are
$\hat {P}=10^6$
Pa,
$(x_c,z_c)=(0,-2000)$
m,
$g=9.81$
m/s
$^2$
,
$c=1450$
m/s,
$h=4000$
m and
$\sigma =200$
m. The error measure is given by
\begin{align} \max _{x,z}\frac {\rho \varPhi _t(x,z,0)+P_0(x,y)}{\hat {P}}, \end{align}
as the strength of the pressure
$\hat {P}$
acts as a linear scaling of the solution.
Tables 2, 3 and 4 present the errors when varying
$N$
,
$k_{\textit{max}}$
and
$\Delta k$
, respectively. When not varied, the default parameter values
$N=100$
,
$k_{\textit{max}}=0.2$
and
$\Delta k=0.0002$
are used. We note that the error given by the default parameters is sufficient for visual convergence. Of the three parameters, the number of included modes
$N$
(table 2) has the most significant effect on the accuracy of the solution because it determines how well the initial condition can be resolved in the vertical direction. The choice of
$k_{\textit{max}}$
(table 3) is affected by the choice of pressure distribution, in our case, leading to exponential decay in the integrand. Finally, all tested discretisation spacings
$\Delta k$
(table 4) resulted in similar error, indicating that we have sufficient resolution of our integrand to numerically perform the integration.
Initial condition reconstruction error for different
$k_{\textit{max}}$
, and
$N=100$
and
$\Delta k=0.0002$
.
Initial condition reconstruction error for different
$\Delta k$
, and
$k_{\textit{max}}=0.2$
and
$N=100$
.
Appendix D. Normalisation coefficients in the case of static compression
In this section, we compute the normalisation coefficients
$D_{\!p}(k)$
needed in the solution in the case of static compression, as discussed in § 5. We begin by introducing a new inner product over the vertical ordinate
\begin{equation} \langle u_1,u_2\rangle _h=\int _{-h}^0 e^{-\gamma z}u_1(z)\overline {u_2(z)}\,\textrm {d}z+\frac {c^2}{g}u_1(0)\overline {u_2(0)}, \end{equation}
and seek to compute
$\langle f_n(\boldsymbol{\cdot };k),f_m(\boldsymbol{\cdot };\tilde {k})\rangle _h$
. Following some tedious calculations, we find
\begin{align} &\left (\mu _n^2-\frac {\gamma ^2}{4}\right )\int _{-h}^0 e^{-\gamma z}f_n(z;k)f_m(z;\tilde {k})\,\textrm {d}z=\big [e^{-\gamma z}f_n'(z;k)f_m(z;\tilde {k})\big ]_{-h}^0-\big [e^{-\gamma z}f_n(z;k)f_m'(z;\tilde {k})\big ]_{-h}^0\nonumber \\ &\quad +\int _{-h}^0 e^{-\gamma z}f_n(z;k)\left (\tilde {\mu }_m^2-\frac {\gamma ^2}{4}\right )f_m(z;\tilde {k})\,\textrm {d}z, \end{align}
which yields
\begin{align} \int _{-h}^0 e^{-\gamma z}f_n(z;k)f_m(z;\tilde {k})\,\textrm {d}z=\frac {\omega _n^2(k)-\omega _m^2(\tilde {k})}{g \left (\mu _n^2-\tilde {\mu }_m^2\right )}. \end{align}
Now, using (D1),
\begin{align} \langle f_n(z;k),f_m(z;\tilde {k})\rangle _h&=\frac {\omega _n^2(k)-\omega _m^2(\tilde {k})}{g \left (\mu _n^2-\tilde {\mu }_m^2\right )}f_n(0;k)f_m(0;\tilde {k})+\frac {c^2}{g},\nonumber \\ &=\frac {c^2}{g}\frac {k^2-\tilde {k}^2}{\mu _n^2-\tilde {\mu }_m^2}. \end{align}
When
$n=m$
, the above inner product reduces to
\begin{align} \lim _{k\rightarrow \tilde {k}} \langle f_m(z;k),f_m(z;\tilde {k})\rangle _h = \frac {c^2}{g}\frac {\tilde {k}}{\tilde {\mu }_m\partial _k\tilde {\mu }_m}, \end{align}
where
\begin{align} & \hspace{1cm}\!\!\! \partial _k\mu _n=\frac {\displaystyle 2\frac {c^2}{g}\left (1-{\gamma }/{2}({\tanh {\mu _n h}}/{\mu _n})\right )k}{\mathcal{A}(\mu _n,k)\tanh {\mu _n h} + \mathcal{B}(\mu _n,k)\mathrm{sech}^2{\mu _n h}+2 {c^2}/{g}\mu _n}, \\[-12pt]\nonumber \end{align}
with
\begin{align} &\hspace{-1cm}\qquad \qquad \mathcal{A}(\mu _n,k)= \left (1-\frac {\gamma c^2}{2g}\right ) +\left ( \frac {\gamma }{2\mu _n}\right )^2-\frac {\gamma c^2}{2g} \left (\left (\frac {k}{\mu _n}\right )^2+\left (\frac {\gamma }{2\mu _n}\right )^2 \right )\!,\nonumber \\ &\hspace{-1cm}\qquad \qquad \mathcal{B}(\mu _n,k) = \mu _nh\left [1- \left ( \frac {\gamma }{2\mu _n}\right )^2 +\frac {\gamma c^2}{2g}\left (\left (\frac {k}{\mu _n}\right )^2+\left (\frac {\gamma }{2\mu _n}\right )^2-1\right ) \right ]\!. \end{align}
Consequently, the inner product
$\langle \phi _n^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k),\phi _{\!p}^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },\tilde {k})\rangle _{\mathcal{E}}$
is computed as
\begin{align} \langle \phi _n^\pm , \phi _{\!p}^\pm \rangle _{\mathcal{E}} &=\int _{-\infty }^\infty \int _{-h}^0 \phi _n^\pm (x,z;k)\overline {\phi _{\!p}^\pm (x,z;\tilde {k})}\,\textrm {d}z\,\textrm {d}x +\frac {c^2}{g}\int _{-\infty }^\infty \phi _n^\pm (x,0;k) \overline {\phi _{\!p}^\pm (x,0;\tilde {k})}\,\textrm {d}x,\nonumber \\ &=\int _{-\infty }^\infty \langle f_n(z;k),f_{\!p}(z;\tilde {k})\rangle _h e^{\mp \textrm {i}(k-\tilde {k}) x}\,\textrm {d}x = 2\pi \frac {c^2}{g}\frac {\tilde {k}}{\tilde {\mu }_n\partial _k\tilde {\mu }_n}\delta (k-\tilde {k})\delta _{np}. \end{align}
Note that in the case of counter-propagating waves, we compute
$\langle \phi _n^\pm (\boldsymbol{\cdot },\boldsymbol{\cdot },k),\phi _{\!p}^\mp (\boldsymbol{\cdot },\boldsymbol{\cdot },\tilde {k})\rangle _{\mathcal{E}}=0$
, which is analogous to (A2f
).
'/%3e%3cpath%20d='M17%207.971L15.711%206.677%2010.912%2011.495%2010.912%200%209.088%200%209.088%2011.495%204.289%206.677%203%207.971%2010%2015z'%20transform='translate(8%208)'/%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
