1. Introduction
Tsunami waves generated from submarine earthquakes or landslides encounter numerous topographies of the ocean floor, such as continental shelves, continental slopes, abyssal planes, ocean trenches and plateaus and mid-oceanic ridges, that lead to the evolution of transoceanic wave propagation (Kulikov et al. Reference Kulikov, Gusiakov, Ivanova and Baranov2016; Salaree & Okal Reference Salaree and Okal2020). The pressure wave generated around the initial disturbance also transforms across different topographies and is entangled with the surface profile. The resulting wave evolution and the underwater pressure wave structure quantify the intensity of the tsunami waves as they make landfall along the coastal water, leading to the severity of impending damages (Miyashita, Mori & Gómez-Ramos Reference Miyashita, Mori and Gómez-Ramos2022). According to the World Health Organization, the estimate of deaths caused by tsunamis between
$1998\;{\rm and}\;2017$
reached
$250\,000$
, with the
$2004$
India Ocean Tsunami alone responsible for more than
$2\,27\,000$
. The World Bank estimated a total loss of more than US
$\$ 7$
billion (http://documents.worldbank.org/curated/en/194061468258308532/World-Bank-response-to-the-Tsunami-disaster). Clearly, the estimation of the tsunami wave profile and pressure wave signature across the ocean, which is influenced by different seafloor topography, is of paramount importance. Here, our objective is to build a mathematical model of submarine earthquake-generated surface waves impacted by arbitrary ocean floor topography in their propagation course, emphasising the computational accuracy of the pressure wave profile that involves contributions from high frequency. Such computational accuracy is lacking in previous calculations, and the present work presents the benchmark computations for pressure wave profiles.
Given the length scale of tsunami wave propagation, the surface tension and viscosity of ocean water are not as important as the ocean compressibility, which paves the way for a mathematical structure associated with the acoustic-gravity wave (AGW) – a low-frequency hydroacoustic wave influenced by gravity, which plays an important role in accurate computations of the pressure inside the water column. These AGW modes generated in a compressible ocean are analogous to the seismic waves in the Earth. They are compressional waveguide modes in the ocean. We note that the first AGW mode becomes the tsunami wave mode when compression is taken to be infinite. For realistic parameters, this mode is similar to the incompressible tsunami wave. Although the mathematical technique to add ocean water compressibility into the physical system is well established through the works of Longuet-Higgins (Reference Longuet-Higgins1950), Miyoshi (Reference Miyoshi1954), Sells (Reference Sells1965) and Yamamoto (Reference Yamamoto1982), a renewed interest in applying the theory to the tsunami wave propagation model is later found through a series of works by Nosov (Nosov Reference Nosov1999; Nosov & Skachko Reference Nosov and Skachko2001; Nosov & Kolesov Reference Nosov and Kolesov2007), which turned out to be influential in the aftermath of the devastating Indian Ocean tsunami. These works pointed out the enhanced accuracy of tsunami wave profile estimation primarily due to the impact of AGW in a compressible ocean compared with the same for an incompressible ocean, a common assumption made in the majority of the works. Later, Stiassnie (Reference Stiassnie2010) and Hendin & Stiassnie (Reference Hendin and Stiassnie2013) provided a solution technique to address the far-field behaviour of the tsunami and AGW modes. With this idea, subsequent researchers have investigated different properties of AGW, including two interacting gravity waves producing AGW(Kadri & Stiassnie Reference Kadri and Stiassnie2013), triad resonances (Kadri Reference Kadri2016; Kadri & Akylas Reference Kadri and Akylas2016), deep-water transport (Kadri Reference Kadri2014), ice shelf breaking (Kadri & Stiassnie Reference Kadri and Stiassnie2012) and wave dissipation on sea ice (Chen, Gilbert & Guyenne Reference Chen, Gilbert and Guyenne2019) or the influence of different ocean floor characteristics such as porous (Chierici, Pignagnoli & Embriaco Reference Chierici, Pignagnoli and Embriaco2010), viscous (Abdolali, Kirby & Bellotti Reference Abdolali, Kirby and Bellotti2015) or rigid but variable bathymetry using depth-integrated mild-slope equations (Sammarco et al. Reference Sammarco, Cecioni, Bellotti and Abdolali2013) or Fourier-series expansion (Das & Meylan Reference Das and Meylan2023b ). However, all of these works neglected one aspect of ocean water compressibility, namely static compression.
Static compression refers to the increase in the density of ocean water as a result of the pressure of the overlying fluid. It is only important for deep ocean water. The associated theory was first developed in the work of Longuet-Higgins (Reference Longuet-Higgins1950) but is often neglected. It has only a small significance for typical ocean depths, and inclusion requires an increase in mathematical difficulty. For example, static compression equations appear in Omira et al. (Reference Omira, Ramalho, Kim, González, Kadri, Miranda, Carrilho and Baptista2022) as part of their analysis of tsunami waves generated by the Hunga Tonga–Hunga Ha‘apai volcanic eruption, but they then neglect this term in their solution. However, a few recent works (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) have shown its importance in the pure gravity mode (tsunami), although their primary focus was on finding the dispersion relation and the propagation speed. In particular, Abdolali & Kirby (Reference Abdolali and Kirby2017) emphasised the static compression-influenced amplification of the phase speed at a higher water depth. Later, Meza-Valle, Kadri & Ortega (Reference Meza-Valle, Kadri and Ortega2023) introduced a technique to obtain the far-field behaviour of waves based on a simplified dispersion relation associated with a large water depth. The primary issue of dealing with static compression at finite depth is the establishment of an inner product between the complicated depth-dependent functions so that more realistic ocean floor topography can be handled. The issue was addressed by Das & Meylan (Reference Das and Meylan2023a ), who provided a complete solution for the canonical problem of vertically moving ocean floor-generated surface wave motion using the eigenfunction expansion method. The simplicity of the method in tackling the mathematical rigour concerning static compression over standard transformation methods makes the approach novel. The method was later extended to three dimensions (Das & Meylan Reference Das and Meylan2024) and to model the transformation of surface waves over a step (Das, Pethiyagoda & Meylan Reference Das, Pethiyagoda and Meylan2024) and a trench (Pethiyagoda et al. Reference Pethiyagoda, Das and Meylan2024b ). Static compression has also been introduced to the propagation of wave propagation models from surface disturbances (Pethiyagoda, Das & Meylan Reference Pethiyagoda, Das and Meylan2025, Reference Pethiyagoda, Das, Bonham and Meylan2024a ), which supplemented the work of Renzi & Dias (Reference Renzi and Dias2014) where static compression was ignored and Liu & Higuera (Reference Liu and Higuera2022) where an incompressible water was assumed.
In reality, the ocean bed shows features of a non-rigid nature, and prior studies emphasised the impact of elastic (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), porous (Gomez & Kadri Reference Gomez and Kadri2023) and viscoelastic (Bahena-Jimenez, Bautista & Méndez Reference Bahena-Jimenez, Bautista and Méndez2023) properties. These features lead to dissipation of the AGW modes, especially the elastic nature that contributes to the vanishing of the leading-order AGW at shallower depth, where it subsequently gets translated to Scholte and Rayleigh waves (Eyov et al. Reference Eyov, Klar, Kadri and Stiassnie2013; Williams & Kadri Reference Williams and Kadri2023) and generation of a very small-amplitude surface-gravity mode (with cutoff frequency) having propagation speed close to the speed of sound (Williams & Kadri Reference Williams and Kadri2023). In addition, due to the temperature gradient and the salinity of the vertical ocean water profile, water density varies in a one-to-one correspondence with the sound speed within. Recent studies have tried to capture the evolution of AGW propagation in such a system (Michele & Renzi Reference Michele and Renzi2020) and also under Lagrangian formulation (Dubois et al. Reference Dubois, Imperiale, Mangeney, Bouchut and Sainte-Marie2023). However, we are focused solely on the topographical properties of the seabed, instead of the nature of it and variable ocean density (consequently, a vertically changing sound speed profile), which requires further rigorous mathematical treatment. We build up the solution of the canonical problem in two dimensions of the rigid seabed, focusing explicitly on the computational accuracy of the pressure wave profile and its transient nature entangled with the associated surface-gravity wave propagation while solving the full dispersion relation, something overlooked in the majority of the prior research, where the far-field solutions are sought.
While the free-surface elevation is dominated by the first AGW mode, the accurate computation of the water-column pressure wave propagation requires higher AGW modes involving high frequencies, which are computationally difficult. The visualisation of AGW through pressure wave propagation inside the ocean water column was not achieved until the recent works of the authors of this paper. The present results calculate this wave in unprecedented detail. Detection of waves gives a great deal of information about the underlying system in many branches of physics, not least earthquakes. Similarly, we believe that detecting these AGWs will give very important information about the associated tsunami wave. Recently, Kadri, Abdolali & Kirby (Reference Kadri, Abdolali and Kirby2025) attempted a similar problem but focused on the quantum tunnelling effect of an AGW, whereas we focus on the effect of the transient initial condition. In this work, we attempt novel computations involving such high-frequency modes for improved pressure wave solutions. Then we are able to extend the solution to tackle arbitrary ocean floor topography using step-type approximations, where the eigenfunction expansion method coupled with the Fourier transformation is used. A few realistic geometries, such as continental slopes and mountain ridges, are taken as test cases. Moreover, a nonlinear temporal growth function associated with ocean floor motion is attempted in addition to a linear one. The manuscript is structured as follows. Section 2 presents the two-dimensional mathematical formulation of the physical problem, followed by the solution method that involves eigenfunction matching across the boundaries between regions in § 3. The time-domain simulations are detailed in § 4. The manuscript concludes with the limitation and future scope in § 5, followed by a brief summary in § 6.
2. Mathematical formulation and modelling assumptions
We briefly discuss the modelling assumptions which underlie this work. The starting point for our work here is the papers by Stiassnie and coworkers (Stiassnie Reference Stiassnie2010; Hendin & Stiassnie Reference Hendin and Stiassnie2013), and almost all of our assumptions follow from those used in this work. We note these same assumptions are also made in other recent work (Kadri et al. Reference Kadri, Abdolali and Kirby2025). We acknowledge that these works did not include the effect of sea floor elasticity, and neither does the work here. We note that this effect may be important (Abdolali et al. Reference Abdolali, Kadri and Kirby2019; Williams & Kadri Reference Williams and Kadri2023) in shallow water. The only modelling change we make to the work of Stiassnie is to include static compression, as formulated in the classical work of Sells (Reference Sells1965). In the literature, there has been discussion about the importance of this term. We include this term primarily to test its effect and because it adds no run time to the numerical solution. Our results here and in previous works show that its effect is small. We note that more complicated models for static compression have been developed (Michele & Renzi Reference Michele and Renzi2020), but we anticipate that they will also only have a small effect. Finally, we add a brief note concerning solving of an equation. We believe that the solution needs to be computed and visualised. We note that a full numerical solution was not presented in Hendin & Stiassnie (Reference Hendin and Stiassnie2013) and many related works. We believe that it is essential to capture the entire pressure throughout the fluid, not simply the motion of the free surface.
Schematic of the physical problem depicting the ocean floor’s vertical movement, which is amplified for visual clarity. The block width and ocean depth have length scales of kilometres, whereas
$a_0$
is in metres.
The generation of free-surface-gravity waves from the finite-time arbitrary vertical movement of a block of the ocean bottom of fixed length in a compressible ocean, including the effect of static background compression, is considered. The vertical movement has a maximum uplift of
$a_0$
and a duration of
$\tau$
seconds, which induces a rise in the free surface and results in wave propagation. The ocean bathymetry is arbitrary and approximated with step-type topography at the locations
$x=l_i\,(i=0,1,\ldots ,m_1)$
with assumptions
$l_0\to -\infty$
and
$l_{m_1}\to \infty$
and local ocean depth
$h_i$
in the region
$(l_{i-1},\,l_i)$
(say region
$\mathcal{R}_i$
). In addition,
$a_0\ll (h_i-h_{\!j})$
(for
$i,j=0,1,\ldots , m_1$
) is imposed to show the significant local ocean depth. Without loss of generality, we assume that the movement takes place in
$\mathcal{R}_k$
(figure 1). A two-dimensional Cartesian coordinate system, centred in the middle of the movement, is considered for a mathematical formulation involving the
$z$
-axis being perpendicular to the horizontal
$x$
-axis and pointing upwards. Under the assumption of irrotational flow, we solve the boundary value problem (BVP) constructed with the velocity potential
$\varPhi (x,z,t)$
inside the fluid domain that takes different forms based on the ocean depth. Assuming a ridge impermeable bottom, the fluid velocity in the vertical direction is given by the bathymetry velocity,
$\mathcal{I}(x,t)$
.
The BVP inside the trench in the time domain is given by
\begin{align} \varPhi ^{(i)}_{xx}+\varPhi ^{(i)}_{zz}&= \frac {1}{c^2}\big (\varPhi ^{(i)}_{tt}+g \varPhi ^{(i)}_{z}\big ),\quad &-h_i\lt z\lt 0,\,x\in \mathcal{R}_i, \\[-12pt]\nonumber \end{align}
\begin{align} \varPhi ^{(i)}_{tt}+g\varPhi ^{(i)}_{z}&=0, \quad &z=0,\,x\in \mathcal{R}_i, \\[-12pt]\nonumber \end{align}
\begin{align} \varPhi ^{(i)}_{z}&=\mathcal{I}(x,t), &z=-h_i,\,x\in \mathcal{R}_i, \end{align}
where
$i=0,\,1,\ldots ,m_1$
,
$c$
represents the speed of sound in the water,
$g$
is the acceleration due to gravity and the vertical velocity along the bottom is assumed separable and given by
\begin{align} \mathcal{I}(x,t)=\left \{\begin{matrix} \mathcal{L}(t) & \mbox{when }x\in \mathcal{R}_k,\\[3pt] 0 & \mbox{when }x\notin \mathcal{R}_k,\end{matrix} \right . \end{align}
where
$\mathcal{L}(t)$
is the temporal velocity of the ocean floor. Note that the term
$g\phi _z^{(i)}$
in ((2.1a
), (2.1b
)) is significant for surface-gravity waves as opposed to the AGW (Michele & Renzi Reference Michele and Renzi2020). However, we are not interested here in the individual impact of either the AGW modes or the surface-gravity mode, but rather the combined effect of both types on the pressure wave propagation, focusing on accurately capturing the qualitative behaviour using higher-order modes. While a decoupling of the BVP based on the mode types (i.e. with or without the terms
$g\phi _z^{(i)}$
) seems the logical choice, we choose to retain the term in the subsequent solution process to maintain a unified set-up.
The matching conditions at the edge of depth change are specified as
\begin{align} \left (\begin{matrix} \varPhi ^{(i)}\\[3pt] \varPhi ^{(i)}_x \end{matrix} \right ) = \left (\begin{matrix} \varPhi ^{(i+1)}\\[3pt] \varPhi ^{(i+1)}_x \end{matrix} \right ) \quad \mbox{at} \quad x=l_i,-h_{i,m}\lt z\lt 0\quad \text{for } \quad m=1,\,2,\ldots ,m_1-1, \end{align}
where
$h_{i,m}=\text{min} (h_i,\,h_{i+1})$
, with the no-flux condition
\begin{align} \left .\begin{matrix} \partial _x\varPhi ^{(i+1)}=0 & \mbox{when }h_i\lt h_{i+1}\\[3pt] \partial _x\varPhi ^{(i)}=0 & \mbox{when }h_i\gt h_{i+1} \end{matrix}\right \}&\mbox{at } x= l_i,\,-h_{i,M}\lt z\lt -h_{i,m}\quad \text{for } \quad i=1,\,2,\ldots ,m_1-1, \end{align}
where
$h_{i,M}=\text{max} (h_i,\,h_{i+1})$
.
3. Solution technique
We apply a Fourier transformation in time of the form
\begin{align} \phi (x,z,\omega )=\int _{-\infty }^{ \infty } \varPhi (x,z,t)e^{- {i}\omega t}{\rm d}t, \end{align}
whose inverse transformation is
\begin{align} \varPhi (x,z,t)=\frac {1}{2\pi }\int _{-\infty }^{\infty } \phi (x,z,\omega )e^{ {i}\omega t}{\rm d}\omega . \end{align}
Note that, in our problem, the solution will be zero for
$t\lt 0$
. This makes the solution here very close to a Laplace transform solution with the variable
$s$
replaced by
$-\mathrm{i}\omega$
, where
$\omega$
is the angular frequency, and the Bromwich inversion contour on the imaginary (
$s$
axis). We note that the Fourier transform is generally more numerically stable than the Laplace transform.
We transform the BVPs, matching conditions and no-flux conditions stated in (2.1)–(2.4) into the frequency-domain BVP in the inner region as
\begin{align} \left (\partial _{xx}+\partial _{zz} \right )\phi ^{(i)}&= \frac {1}{c^2}\big (-\omega ^2+g\partial _{z}\big )\phi ^{(i) },\quad &h_i\lt z\lt 0,\,x\in \mathcal{R}_i, \\[-12pt]\nonumber \end{align}
\begin{align} \big (-\omega ^2+g\partial _{z}\big ) \phi ^{(i)}&= 0,\quad &z=0,\,x\in \mathcal{R}_i, \\[-12pt]\nonumber \end{align}
\begin{align} \partial _z\phi ^{(i)}&=\xi (\omega ), \quad &z=-h_i,\,x\in \mathcal{R}_i, \end{align}
where
\begin{align} \xi (\omega )=\left \{\begin{matrix} \displaystyle \int _{-\infty }^{ \infty } \mathcal{L}(t)e^{-{i}\omega t}{\rm d}t & \mbox{when }x\in \mathcal{R}_k,\\[3pt] 0 & \mbox{when }x\notin \mathcal{R}_k. \end{matrix} \right . \end{align}
The matching conditions at
$x=l_i$
are transformed into
\begin{align} \left (\begin{matrix} \phi ^{(i)}\\[3pt] \phi ^{(i)}_x \end{matrix} \right ) = \left (\begin{matrix} \phi ^{(i+1)}\\ \phi ^{(i+1)}_x \end{matrix} \right ) &\quad \mbox{along} \quad h_{i,m}\lt z\lt 0, \\[-12pt]\nonumber \end{align}
\begin{align} \left .\begin{matrix} \partial _x\phi ^{(i+1)}=0 & \mbox{when }h_i\lt h_{i+1}\\[3pt] \partial _x\phi ^{(i)}=0 & \mbox{when }h_i\gt h_{i+1} \end{matrix}\right \}&\quad \mbox{along}\quad -h_{i,M}\lt z\lt -h_{i,m}. \end{align}
Letting
$\varphi (x,z,\omega )$
be the corresponding potential function such that
\begin{align} \phi ^{(i)}(x,z,\omega )=\xi (\omega ) \varphi ^{(i)}(x,z,\omega ), \end{align}
we can write the potential as
\begin{align} \varphi ^{(i)}&=\sum _{n=0}^{\infty } \left \{b_{i,n}e^{{i}k_n^{(i)} (x-l_{i-1})} + c_{i,n}e^{-{i}k_n^{(i)} (x-l_i)} \right \}f_n^{(i)}(z)+f^{p(i)}(z),\quad \mbox{for}\quad i=1,\ldots ,m_1, \end{align}
with
$b_{1,n}=0$
,
$c_{m_1,n}=0$
for all
$n$
. The depth-dependent functions are expressed as (following Das & Meylan Reference Das and Meylan2023a
)
\begin{align} f_n^{(i)}(z)=&\frac {\displaystyle e^{(\gamma z/2)}\left ( \frac {\gamma }{2\mu ^{(i)}_n}\sinh {\mu ^{(i)}_n(z+h_i)}-\cosh {\mu ^{(i)}_n(z+h_i)} \right )}{\frac {\gamma }{2\mu ^{(i)}_n}\sinh {( \mu ^{(i)}_n h_i)}-\cosh {(\mu ^{(i)}_n h_i)}}\,\mbox{for } \ i=1,2,\ldots ,m_1, \\[-12pt]\nonumber \end{align}
\begin{align} f^{p(i)}(z)=&\left \{\begin{matrix} \displaystyle \frac {\displaystyle e^{\frac {\gamma (z+h_1)}{2}}\left [\left ( \frac {\omega ^2}{gk_s}-\frac {\gamma }{2k_s} \right )\sinh {(k_sz)} +\cosh {(k_sz)} \right ]}{\displaystyle \left ( \frac {\omega ^2}{c^2k_s} - \frac {\omega ^2\gamma }{2gk_s} \right ) \sinh {(k_sh_1)}+\frac {\omega ^2}{g}\cosh {(k_sh_i)}} & \mbox{for }i=k,\\ 0 & \mbox{for } \ i\neq k, \end{matrix} \right . \end{align}
with
$\gamma =g/c^2$
and
$\mu ^{(i)}_n$
satisfying the following dispersion relation:
\begin{align} \frac {\omega ^2}{g}=\mu ^{(i)}_n \frac {\displaystyle \left [1- \left ( \frac {\gamma }{2\mu ^{(i)}_n} \right )^2 \right ]\tanh {\mu ^{(i)}_n h_i}}{\displaystyle 1-\left (\frac { \gamma }{2\mu ^{(i)}_n} \right )\tanh {\mu ^{(i)}_n h_i}}\quad \mbox{for }i=1,2,\ldots ,m_1, \end{align}
and
\begin{equation} k_s^2=\frac {\omega ^2}{c^2}-\left (\frac {\gamma }{2}\right )^2\quad \mbox{and} \quad (k^{(i)}_n)^2=\frac {\omega ^2}{c^2} -\left (\frac {\gamma }{2}\right )^2+(\mu ^{(i)}_n)^2. \end{equation}
Applying matching conditions at
$x=l_i$
$(i=1,2,\ldots , m_1)$
yields
\begin{align} &\sum _{n=0}^{\infty } \left \{b_{i,n}e^{{i}k_n^{(i)} (l_i-l_{i-1})} + c_{i,n} \right \}f_n^{(i)}(z)+f^{p(i)}(z) = \sum _{n=0}^{\infty } \left \{b_{i+1,n}+ c_{i+1,n}e^{-{i}k_n^{(i+1)} (l_i-l_{i+1})} \right \}\nonumber \\ &\hspace {3.2in}\times f_n^{(i+1)}(z)+f^{p(i+1)}(z), \\[-12pt]\nonumber \end{align}
\begin{align} &\sum _{n=0}^{\infty } k_n^{(i)}\left \{b_{i,n}e^{{i}k_n^{(i)} (l_i-l_{i-1})} - c_{i,n} \right \}f_n^{(i)}(z) = \sum _{n=0}^{\infty }k_n^{(i+1)} \left \{b_{i+1,n}- c_{i+1,n}e^{-{i}k_n^{(i+1)} (l_i-l_{i+1})} \right \}\nonumber \\ &\hspace {2.8in}\times f_n^{(i+1)}(z) .\end{align}
Now, we construct the system of equations.
3.1. Construction of the system of equations
To construct a system of equations to obtain the coefficients
$b_{i,n}$
, and
$c_{i,n}$
after truncating the series up to
$n=N$
, we take the inner product defined as
\begin{align} \left \langle \mathcal{F}_i(z),\, \mathcal{F}_{i+1}(z) \right \rangle _{c}^{h}:= \int _{-h}^0{e^{-\gamma z}\mathcal{F}_i(z)\mathcal{F}_{i+1}(z) \,{\rm d}z}\quad \mbox{for}\quad i=1,2,\ldots ,m_1, \end{align}
of (3.11) with
$\mathcal{F}_i=f_n^{(i)}$
and
$h=h_{i,m}$
or
$h_{i,M}$
. Specifically, we use
$\lt ,\gt _{c}^{h_{i,m}}$
, and
$\lt ,\gt _{c}^{h_{i,M}}$
respectively in ((3.11a
)–(3.11b
)). Using the orthogonality of either
$f_n^{(i)}(z)$
or
$f_n^{(i+1)}(z)$
, we get the following system of equations for each
$j=0,1,\ldots ,N$
:
\begin{align} \sum _{n=0}^{N}&\left \{e^{{i}k_n^{(i)}(l_i-l_{i-1})} b_{i,n}+ c_{i,n}\right \}\left \langle f^{(i)}_n,\, f^{(i+1)}_{\!j} \right \rangle _{c}^{h_{i,m}}+ \left \langle f^{p(i)}, \,f^{(i+1)}_{\!j}\right \rangle _{c }^{h_{i,m}}\nonumber \\ &=\left \{b_{i+1,j}+c_{i+1,j}e^{-{i}k_{\!j}^{(i+1)}(l_i-l_{i+1})} \right \}\left \langle f^{(i+1)}_{\!j},\, f^{(i+1)}_{\!j}\right \rangle _{c }^{h_{i,m}}+ \left \langle f^{p(i+1)}, \,f^{(i+1)}_{\!j}\right \rangle _{c }^{h_{i,m}}, \\[-12pt]\nonumber \end{align}
\begin{align} k_{\!j}^{(i)}&\left \{e^{{i}k_{\!j}^{(i)}(l_i-l_{i-1})} b_{i,j}- c_{i,j} \right \}\left \langle f^{(i)}_{\!j},\, f^{(i)}_{\!j} \right \rangle _{c}^{h_{i,M}}\nonumber \\ &=\sum _{n=0}^{N} k_n^{(i+1)}\left \{b_{i+1,n}-c_{i+1,n}e^{-{i}k_n^{(i+1)}(l_i-l_{i+1})} \right \}\left \langle f^{(i+1)}_n ,\, f^{(i)}_{\!j}\right \rangle _{c }^{h_{i,M}}, \end{align}
when
$h_{i+1}\lt h_i$
, or
\begin{align} &\left \{e^{{i}k_{\!j}^{(i)}(l_i-l_{i-1})} b_{i,j}+ c_{i,j}\right \}\left \langle f^{(i)}_{\!j},\, f^{(i)}_{\!j} \right \rangle _{c}^{h_{i,m}}+ \left \langle f^{p(i)}, \,f^{(i)}_{\!j}\right \rangle _{c }^{h_{i,m}} \nonumber \\ &\qquad =\,\sum _{n=0}^{N}\left \{b_{i+1,n}+c_{i+1,n}e^{-{i}k_n^{(i+1)}(l_i-l_{i+1})} \right \}\left \langle f^{(i+1)}_n,\, f^{(i)}_{\!j}\right \rangle _{c }^{h_{i,m}}+ \left \langle f^{p(i+1)}, \,f^{(i)}_{\!j}\right \rangle _{c }^{h_{i,m}}, \\[-12pt]\nonumber \end{align}
\begin{align} \sum _{n=0}^{N}&k_n^{(i)} \left \{e^{{i}k_n^{(i)}(l_i-l_{i-1})} b_{i,n}- c_{i,n}\right \}\left \langle f^{(i)}_n,\, f^{(i+1)}_{\!j} \right \rangle _{c}^{h_{i,M}} \nonumber \\ &\qquad =\,k_{\!j}^{(i+1)}\left \{b_{i+1,j}-c_{i+1,j}e^{-{i}k_{\!j}^{(i+1)}(l_i-l_{i+1})} \right \}\left \langle f^{(i+1)}_{\!j},\, f^{(i+1)}_{\!j}\right \rangle _{c }^{h_{i,M}}, \end{align}
when
$h_{i+1}\gt h_i$
, where the inner products are given by
\begin{align} \left \langle f^{(i)}_n,\,f^{(i)}_{\!j} \right \rangle _c^{h_i}=&\frac {\displaystyle \frac { \sin {2\mu ^{(i)}_nh_i}}{4\mu ^{(i)}_n}\left (1- \frac {\gamma ^2}{4(\mu ^{(i)}_n)^2} \right )+\frac {h_i}{2}\left (1+\frac { \gamma ^2}{4(\mu ^{(i)}_n)^2}\right )+\frac {\gamma \left ( \cos {2\mu ^{(i)}_nh_i} -1\right )}{4(\mu ^{ (i)}_n)^2}}{\left [\frac {\gamma }{2\mu ^{(i)}_n} \sin {(\mu ^{(i)}_n h_i)}-\cos {(\mu ^{ (i)}_n h_i)} \right ]^2} \delta _{jn}, \end{align}
\begin{equation} \left \langle f^{(1)}_{\!j},\,f^p \right \rangle _c^{h_1}= \frac {f^{(1)}_{\!j}(-h_1)e^{\gamma h_1}}{{\mu ^{(1)}_{\!j}}^2-k_s^2}, \end{equation}
\begin{align} &\hspace{-14.5pc} \left \langle f^{(i)}_{\!j},\,f^{(i+1)}_n\right \rangle _{c}^{h}= \frac {f^{(i)}_{\!j}(-h){f^{(i+1)}_n}'(-h)e^{\gamma h}}{{\mu ^{(i)}_{\!j}}^2-{\mu ^{(i+1)}_n}^2}, \end{align}
where
$\delta _{mn}$
is the Kronecker delta. After scaling by the repeated inner product,
$\langle f^{(i)}_{\!j},\, f^{(i)}_{\!j} \rangle _{c}^{h}$
or
$\langle f^{(i+1)}_{\!j},\, f^{(i+1)}_{\!j} \rangle _{c}^{h}$
, we recover the linear subproblem
\begin{equation} \left [\begin{matrix} N_{i1} & N_{i2} & N_{i3} & N_{i4}\\[3pt] N_{i5} & N_{i6} & N_{i7} & N_{i8} \end{matrix}\right ] \left [\begin{matrix} B_{i} \\[3pt] C_{i}\\[3pt] B_{i+1} \\[3pt] C_{i+1} \end{matrix}\right ] = \left [\begin{matrix} \hat {d}_i\\[3pt] 0_d \end{matrix}\right ]\!, \end{equation}
where the
$(j,n)$
th element of each
$N_{i\kappa }$
is defined as
\begin{align} N_{i1}&=\left \{ \begin{matrix} \displaystyle \left [e^{{i}k_n^{(i)}(l_i-l_{i-1})}\frac {\left \langle f^{(i)}_n,\, f^{(i+1)}_{\!j} \right \rangle _{c}^{h_{i,m}}}{\left \langle f^{(i+1)}_{\!j},\, f^{(i+1)}_{\!j}\right \rangle _{c }^{h_{i,m}}}\right ]_{jn} & \mbox{when } h_{i+1}\lt h_i,\\[17pt] \mbox{diag} \Big (e^{ {i} k_{\!j}^{(1)}(l_i-l_{i-1})} \Big)_{1\times (N+1)} & \mbox{when } h_{i+1}\gt h_i, \end{matrix} \right .\nonumber \\[17pt] N_{i2}&=\left \{ \begin{matrix} \displaystyle \left [\frac {\left \langle f^{(i)}_n,\, f^{(i+1)}_{\!j} \right \rangle _{c}^{h_{i,m}}}{\left \langle f^{(i+1)}_{\!j},\, f^{(i+1)}_{\!j}\right \rangle _{c }^{h_{i,m}}}\right ]_{jn} & \mbox{when } h_{i+1}\lt h_i,\\[15pt] \mbox{diag}\left \{(1)_{1 \times (N+1)}\right \} & \mbox{when } h_{i+1}\gt h_i, \end{matrix} \right .\nonumber \\[17pt] N_{i3}&=\left \{ \begin{matrix} -\mbox{diag}\left \{(1)_{1 \times (N+1)}\right \} & \mbox{when } h_{i+1}\lt h_i,\\[17pt] \displaystyle -\left [\frac {\left \langle f^{(i+1)}_n,\, f^{(i)}_{\!j} \right \rangle _{c}^{h_{i,m}}}{\left \langle f^{(i)}_{\!j},\, f^{(i)}_{\!j}\right \rangle _{c }^{h_{i,m}}}\right ]_{jn} & \mbox{when } h_{i+1}\gt h_i, \end{matrix} \right .\nonumber \\[17pt] N_{i4}&=\left \{ \begin{matrix} -\mbox{diag} \left (e^{-{i} k_{\!j}^{(i+1)}(l_i-l_{i+1})} \right )_{1\times (N+1)} & \mbox{when } h_{i+1}\lt h_i,\\[17pt] -\displaystyle \left [e^{-{i}k_n^{(i+1)}(l_i-l_{i+1})} \frac {\left \langle f^{(i+1)}_n,\, f^{(i)}_{\!j} \right \rangle _{c}^{h_{i,m}}}{\left \langle f^{(i)}_{\!j},\, f^{(i)}_{\!j}\right \rangle _{c }^{h_{i,m}}}\right ]_{jn} & \mbox{when } h_{i+1}\gt h_i, \end{matrix} \right .\nonumber \\[17pt] N_{i5}&=\left \{ \begin{matrix} \mbox{diag} \left (e^{ {i} k_{\!j}^{(i)}(l_i-l_{i-1})} \right )_{1\times (N+1)} & \mbox{when } h_{i+1}\lt h_i,\\[17pt] \displaystyle \left [\frac { k_n^{(i)}e^{{i}k_n^{(i)}(l_i-l_{i-1})}}{k_{\!j}^{(i+1)}} \frac {\left \langle f^{(i)}_n,\, f^{(i+1)}_{\!j} \right \rangle _{c}^{h_{i,M}}}{\left \langle f^{(i+1)}_{\!j},\, f^{(i+1)}_{\!j}\right \rangle _{c }^{h_{i,M}}}\right ]_{jn} & \mbox{when } h_{i+1}\gt h_i, \end{matrix} \right .\nonumber \\[17pt] N_{i6}&=\left \{ \begin{matrix} -\mbox{diag}\left \{(1)_{1 \times (N+1)}\right \} & \mbox{when } h_{i+1}\lt h_i,\\[17pt] \displaystyle -\left [\frac { k_n^{(i)}}{k_{\!j}^{(i+1)}} \frac {\left \langle f^{(i)}_n,\, f^{(i+1)}_{\!j} \right \rangle _{c}^{h_{i,M}}}{\left \langle f^{(i+1)}_{\!j},\, f^{(i+1)}_{\!j}\right \rangle _{c }^{h_{i,M}}}\right ]_{jn} & \mbox{when } h_{i+1}\gt h_i, \end{matrix} \right .\nonumber \\[17pt] N_{i7}&=\left \{ \begin{matrix} -\displaystyle \left [\frac { k_n^{(i+1)}}{k_{\!j}^{(i)}} \frac {\left \langle f^{(i+1)}_n,\, f^{(i)}_{\!j} \right \rangle _{c}^{h_{i,M}}}{\left \langle f^{(i)}_{\!j},\, f^{(i)}_{\!j} \right \rangle _{c }^{h_{i,M}}} \right ]_{jn} & \mbox{when } h_{i+1}\lt h_i,\\[17pt] -\mbox{diag}\left \{(1)_{1 \times (N+1)}\right \} & \mbox{when } h_{i+1}\gt h_i, \end{matrix} \right .\nonumber \\[17pt] N_{i8}&=\left \{ \begin{matrix} \displaystyle \left [\frac { k_n^{(i+1)}e^{-{i}k_n^{(i+1)}(l_i-l_{i+1})}}{k_{\!j}^{(i)}} \frac {\left \langle f^{(i+1)}_n,\, f^{(i)}_{\!j} \right \rangle _{c}^{h_{i,M}}}{\left \langle f^{(i)}_{\!j},\, f^{(i)}_{\!j} \right \rangle _{c }^{h_{i,M}}} \right ]_{jn} & \mbox{when } h_{i+1}\lt h_i,\\[17pt] \mbox{diag} \left (e^{-{i} k_{\!j}^{(i+1)}(l_i-l_{i+1})} \right )_{1\times (N+1)} & \mbox{when } h_{i+1}\gt h_i, \end{matrix} \right . \end{align}
for
$j,n=0,1,\ldots ,N$
, and
$B_{i}=[b_{i,0},b_{i,1},\ldots , b_{i,N}]^T$
,
$C_{i}=[c_{i,0},c_{i,1},\ldots ,c_{i,N}]^T$
are column vectors of length
$N+1$
containing unknowns and
$\hat {d}_{i}=[d_{i,0},d_{i,1},\ldots ,d_{i,N}]^T$
is a known column vector of length
$N+1$
defined by
\begin{align} \hat {d}_{i,j}=\begin{cases} \displaystyle \frac {\left \langle f^{p(i+1)}, \, f^{(i+1)}_{\!j}\right \rangle _{c }^{h_{i,m}}-\left \langle f^{p(i)}, \, f^{(i+1)}_{\!j}\right \rangle _{c }^{h_{i,m}}}{\left \langle f^{(i+1)}_{\!j},\, f^{(i+1)}_{\!j} \right \rangle _{c }^{h_{i,m}}}&h_{i+1}\lt h_i\\[17pt] \displaystyle \frac {\left \langle f^{p(i+1)}, \,f^{(i)}_{\!j} \right \rangle _{c }^{h_{i,m}}-\left \langle f^{p(i)}, \,f^{(i)}_{\!j} \right \rangle _{c }^{h_{i,m}}}{\left \langle f^{(i)}_{\!j},\, f^{(i)}_{\!j} \right \rangle _{c }^{h_{i,m}}}&h_{i+1}\gt h_i ,\end{cases} \end{align}
for
$j=0,1,\ldots ,N$
, and
$0_d$
is a column vector of zeros of length
$N+1$
. We combine the subproblem for every interface into a single linear system, recalling that
$b_{1,n}=0$
,
$c_{m_1,n}=0$
for all
$n$
\begin{align} &\left [\begin{array}{cccccccccccc} N_{1\,2} & N_{1\,3} & N_{1\,4} & {\cdots} & 0_A & 0_A & 0_A & 0_A & {\cdots} & 0_A & 0_A & 0_A\\[11pt] N_{1\,6} & N_{1\,7} & N_{1\,8} & {\cdots} & 0_A & 0_A & 0_A & 0_A & {\cdots} & 0_A & 0_A & 0_A\\[11pt] \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\[11pt] 0_A & 0_A & 0_A & {\cdots} & N_{i1} & N_{i2} & N_{i3} & N_{i4} & {\cdots} & 0_A & 0_A & 0_A\\[11pt] 0_A & 0_A & 0_A & {\cdots} & N_{i5} & N_{i6} & N_{i7} & N_{i8} & {\cdots} & 0_A & 0_A & 0_A\\[11pt] \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\[11pt] 0_A & 0_A & 0_A & {\cdots} & 0_A & 0_A & 0_A & 0_A & {\cdots} & N_{m_1-1\,1} & N_{m_1-1\,2} & N_{m_1-1\,3} \\[11pt] 0_A & 0_A & 0_A & {\cdots} & 0_A & 0_A & 0_A & 0_A & {\cdots} & N_{m_1-1\,5} & N_{m_1-1\,6} & N_{m_1-1\,7} \end{array}\right ]\nonumber \\[7pt] &\qquad \times \left [ \begin{array}{cccccccccccc} C_{1} & B_{2} & C_{2} & {\cdots} & B_{i} & C_{i} & B_{i+1} & C_{i+1} & {\cdots} & B_{m_1-1} & C_{m_1-1} & B_{m_1} \end{array}\right ]^{T} \nonumber \\[7pt] &\qquad =\left [\begin{array}{cccccccccccc} \hat {d}_1 & 0_d & {\cdots} & \hat {d}_i & 0_d & {\cdots} & \hat {d}_{m_1-1} & 0_d \end{array}\right ]^{T}, \end{align}
where each
$0_A$
is a zero square matrix of order
$N+1$
. We solve (3.18) to calculate the unknown coefficient and evaluate the potential functions in (3.7).
Once we get the coefficients, we use (3.6) and then apply the inverse Fourier transformation to obtain the potential functions and, hence, the surface profile and the pressure distribution inside the water column in the time domain.
4. Time-domain simulation
The potential functions in the time domain are obtained by combining (3.1b ), (3.6) and (3.7) as
\begin{align} \varPhi ^{(i)}(x,z,t)=\frac {1}{2 \pi } \int _{-\infty }^{\infty }\xi (\omega ) &\left ( \sum _{n=0}^{\infty }\{ b_{i,n} e^{{i}k_n^{(i)} (|x|-l_{i-1})} + c_{i,n} e^{-{i}k_n^{(i)} (|x|-l_i)}\}f_n^{(i)}(z) + f^{p(i)}(z) \right )\nonumber \\ &\times e^{ {i}\omega t}{\rm d}\omega . \end{align}
The free-surface elevation is computed as
\begin{equation} \eta ^{(i)}={\mathrm{Re}}\left (-\frac {1}{g}\left .\frac { \partial \varPhi ^{(i)}}{\partial t}\right |_{z=0}\right )={\mathrm{Re}}\left (-\frac {{\rm i}}{2 \pi g}\int _{-\infty }^{\infty }\omega \phi ^{(i)}(x,0,\omega )e^{ {i}\omega t}\,{\rm d}\omega \right ) \end{equation}
for
$(i)=1,2,\ldots ,N$
, where
${\mathrm{Re}}(\boldsymbol{\cdot })$
represents the real part of a complex quantity. The pressure distribution inside the water column is computed as
\begin{equation} P^{(i)}=-\rho \frac { \partial \varPhi ^{(i)}}{\partial t}=-\frac {\textrm {i}\rho }{2 \pi }\int _{-\infty }^{\infty } \omega \phi ^{(i)}(x,z,\omega )e^{ {i}\omega t}\,{\rm d}\omega. \end{equation}
All integrals in (4.1)–(4.3) have simple poles along the real axis corresponding to when the determinant of the coefficient matrix in (3.18) is zero. The contour of integration passes over the poles for
$\omega \lt 0$
and under for
$\omega \gt 0$
. The integral is evaluated following Pethiyagoda et al. (Reference Pethiyagoda, Das and Meylan2024b
), whereby the integral over the truncated domain
$\omega \in [0,\omega _{{max}}]$
is subdivided in regions near the poles and elsewhere. The regions near the poles are further modified by adding and subtracting the singularity resulting in continuous integral and a closed form evaluation of a Cauchy principal valued integral and half-residue contribution. All remaining integrals are evaluated numerically using a Filon-type quadrature.
Now, we consider a few specific bathymetry profiles and two different velocity profiles for bottom rupture, constant velocity
$\mathcal{L}(t)=\displaystyle ({a_0}/{\tau })H(\tau -t)H(t)$
and parabolic velocity
$\mathcal{L}(t)=\displaystyle({6a_0}/{\tau })t(\tau -t)H(\tau -t)H(t)$
, where
$H(t)$
is the Heaviside step function. The fixed parameter values are
$a_0=1$
m,
$c=1450$
m s−1,
$g=9.8$
m s−
$^2$
,
$\tau =10$
s and
$N=30$
, unless otherwise mentioned. Our choice of the topography parameters in the following subsections is idealistic, and is utilised to show the robustness of the method. More accurate topography from the observational data, something we have not collected at present, may be fit using a proper step approximation to validate the solution with observational data. Nevertheless, we move forward with the idealistic parameters. For the calculations here, there are no solutions in the literature to validate against. However, the solution is mathematically equivalent to our previous solutions for the constant depth case (Das & Meylan Reference Das and Meylan2024), which were validated against a Green’s function formulation. We have confidence in our results because we can check that the frequency-domain solutions satisfy the boundary conditions (noting they must satisfy the Helmholtz equation by construction).
4.1. Continental shelf and slope
In this subsection, we show the surface profile and pressure distribution inside the water column for waves generated near continental shelf- and slope-type bathymetries. Our continental shelf is characterised by a maximum depth of
$3.5$
km and a minimum depth of
$500$
m. The change in depth occurs over five equally spaced steps from
$x=40$
to
$x=100$
km (figure 2), leading to an average inclination of approximately
$2.3^\circ$
. The region
$x\in [0,20]$
km is chosen as our rupture zone.
4.1.1. Constant rupture velocity
Compression waves near a continental shelf caused by a constant velocity rupture (Movie 1 shows full animation are available at https:/doi,org/10.1017/jfm.2026.11821). Each time panel presents the free-surface profile above a heat map of the dynamic pressure in the fluid. The red dotted lines in the free-surface plot indicate the edge of the rupture zone. The rupture surface is indicated in the pressure plot by the red line from
$0\lt x\lt 20$
.
Figure 2. Long description
Panel A: A heatmap plot displays pressure wave propagation at 2.25 seconds. The x-axis ranges from -150 to 150 kilometres, and the z-axis ranges from 0 to -3.5 kilometres. The colour scale indicates pressure in Pascals, with blue representing lower pressure and yellow representing higher pressure. The surface wave profile above the heatmap plot shows variations in wave height over the same x-axis range. Panel B and C: Equivalent heatmap plot and surface wave profiles show pressure wave propagation and variations in wave height at 10.25 and 78.25 seconds, respectively. The axes and colour scale are the same as in Panel A.
Figure 3. Long description
Panel A: A heatmap plot displays pressure wave propagation at 165.25 seconds. The x-axis ranges from -150 to 150 kilometres, and the z-axis ranges from 0 to -3.5 kilometres. The colour scale indicates pressure in Pascals, with blue representing lower pressure and yellow representing higher pressure. The surface wave profile above the heatmap plot shows variations in wave height over the same x-axis range. Panel B: Another heatmap plot and surface wave profile shows pressure wave propagation and variations in wave height at 675.25 seconds. The axes and colour scale are the same as in Panel A.
For a constant velocity rupture, the evolution of the surface wave profile is shown in figures 2 and 3 with snapshots at
$t=2.25,\,10.25,\,78.25,\,165.25$
and
$675.25$
s. In the first snapshot, we see a column of positive pressure rising from the rupture surface, travelling at the speed of sound, driven by the constant-speed movement of the boundary. This pressure column will reflect off of the water’s surface, changing its sign before continuing back down to the seabed and reflecting again without changing its sign. As this process repeats, the pressure column makes four trips between the seabed and the surface such that it is positive and moving up when the rupture zone stops moving. The combined interference of this pressure column results in a localised positive pressure that will continue to bounce between the seabed and the surface (figure 2
b).
As time continues, pressure waves peel off from the localised pressure pulse and propagate away from the rupture region (figure 2 c). These pressure waves slowly align themselves with the vertical axis as they propagate away from the continental shelf, in the constant depth region, because the finite-depth domain acts as a waveguide. The pressure waves travelling towards the continental shelf scatter at the depth changes and result in a higher-frequency pressure pattern over the shallower depth regions (figure 3).
Focussing on the free-surface profiles, we can see that the effect of compressibility means that the surface only begins to move noticeably after the pressure column has had time to reach the surface (figure 2). Once the surface starts to move, it lifts above the rupture (figure 2 b) before splitting into two gravity waves moving in both directions, which is consistent with early-time symmetric solutions described in the author’s previous works (Das et al. Reference Das, Pethiyagoda and Meylan2024; Pethiyagoda et al. Reference Pethiyagoda, Das and Meylan2024b ). Comparing the two gravity waves at a later time point (figure 3), we can see that the decrease in depth on the right has caused the travelling wave to become steeper. Additionally, the previously mentioned pressure waves generated the AGWs, which appear as small ripples on the surface. Movie 1 and Movie 2 show the simulation videos for shorter and longer time scales.
Compression waves near a continental shelf caused by a parabolic rupture velocity (check Movie 3 for the full animation). Each time panel presents the free-surface profile above a heat map of the dynamic pressure in the fluid. The red dotted lines in the free-surface plot indicate the edge of the rupture zone. The rupture surface is indicated in the pressure plot by the red line from
$0\lt x\lt 20$
.
Figure 4. Long description
Panel A: A heatmap plot displays pressure wave propagation at 10.25 seconds. The x-axis ranges from -150 to 150 kilometres, and the z-axis ranges from 0 to -3.5 kilometres. The colour scale indicates pressure in Pascals, with blue representing lower pressure and yellow representing higher pressure. The surface wave profile above the heatmap plot shows variations in wave height over the same x-axis range. Panel B and C: Equivalent heatmap plot and surface wave profiles show pressure wave propagation and variations in wave height at 165.25 and 675.25 seconds, respectively. The axes and colour scale are the same as in Panel A.
4.1.2. Parabolic rupture velocity
The solutions for a parabolic rupture velocity are shown in figure 4 with snapshots at
$t=10.25,\,165.25$
and
$675.25$
s. The dynamic pressure above the rupture exhibits a more non-localised profile in the vertical direction compared with the previous constant velocity rupture solutions (figure 3). This profile leads to a more defined banded structure in the pressure waves propagating to the left away from the rupture (figure 4
b). As a consequence, the AGWs on the surface are of a greater magnitude than was previously seen for the constant velocity rupture.
The initial banded structure (i.e. low frequency in the
$z$
-direction) of the right-propagating pressure waves leads to greater reflection off of the continental shelf and, therefore, significantly smaller AGWs appearing on the continental shelf when compared with the constant velocity rupture solutions. Additionally, the reflected pressure waves interact with original left-propagating waves, leading to much larger AGWs in both wavelength and amplitude compared with the early-time solutions (figure 4
c). Movie 3 shows the simulation video.
4.2. Trench and ridge
In this subsection, we show the surface profile and dynamic pressure inside the water column for waves generated near an underwater trench mountain ridge. The bathymetry in kilometres is taken as a piecewise constant approximation of
\begin{align} h(x)=-2.5-6(x-20)(x-80)H(x-20)H(x-80)+6(x-120)(x-180)H(x-120)H(x-180). \end{align}
The parabolic arcs that define the trench and ridge are divided into five levels (figure 5). The region
$x\in [-10,10]$
km is chosen as our rupture zone.
4.2.1. Constant rupture velocity
Compression waves near an underwater trench and a mountain ridge caused by a constant velocity rupture. Each time panel presents the free-surface profile above a heat map of the dynamic pressure in the fluid. The red dotted lines in the free-surface plot indicate the edge of the rupture zone. The rupture surface is indicated in the pressure plot by the red line from
$-10\lt x\lt 10$
. Movie 4 (a) and Movie 4 (b) show full animations.
Figure 5. Long description
Panel A: A heatmap plot displays pressure wave propagation at 54 seconds. The x-axis ranges from -40 to 200 kilometres, and the z-axis ranges from 0 to -4 kilometres. The colour scale indicates pressure in Pascals, with blue representing lower pressure and yellow representing higher pressure. The surface wave profile above the heatmap plot shows variations in wave height over the same x-axis range. Panel B and C: Equivalent heatmap plot and surface wave profiles show pressure wave propagation and variations in wave height at 312 and 943 seconds, respectively. The axes and colour scale are the same as in Panel A.
Figure 5 provides snapshots of the surface displacement and dynamic pressure caused by a constant velocity rupture for
$t=54,\,312$
and
$943$
s. The first thing to address is that, even though the rupture speed is the same as that presented in § 4.1.1, the shallower rupture location has resulted in different interference patterns in the pressure and, therefore, a less localised pressure pulse bouncing between the surface and the bottom. This change in local pressure behaviour has the effect of altering the pressure waves that propagate away from the rupture region, leading them to appear more banded than those observed in figure 2 but still less so than those created by the parabolic rupture velocity (figure 4).
As with previous examples, pressure waves begin propagating away from the rupture location. These waves move over the trench with some reflection occurring as the waves leave the trench. These reflections begin to excite a standing wave pattern within the trench. As the waves move past the trench, they interact with and mostly reflect back off of the ridge, much like the waves in § 4.1.2. These additional reflected waves from the ridge return past the trench and further excite the standing waves within the trench such that they are readily apparent on the free surface (figure 5 c).
The gravity waves generated by the rupture behave in much the same way as the previous examples. The gravity waves propagate away from the rupture area, maintaining the same qualitative shape while slightly broadening as they pass over the trench and steepening as they pass over the ridge. The simulations are provided as Movie 4 (a) and Movie 4 (b).
4.2.2. Parabolic rupture velocity
Compression waves near an underwater trench and a mountain ridge caused by a parabolic rupture velocity. Each time panel presents the free-surface profile above a heat map of the dynamic pressure in the fluid. The red dotted lines in the free-surface plot indicate the edge of the rupture zone. The rupture surface is indicated in the pressure plot by the red line from
$-10\lt x\lt 10$
. Movie 5 (a) and Movie 5 (b) show full animations.
Figure 6. Long description
Panel A: A heatmap plot displays pressure wave propagation at 54 seconds. The x-axis ranges from -40 to 200 kilometres, and the z-axis ranges from 0 to -4 kilometres. The colour scale indicates pressure in Pascals, with blue representing lower pressure and yellow representing higher pressure. The surface wave profile above the heatmap plot shows variations in wave height over the same x-axis range. Panel B and C: Equivalent heatmap plot and surface wave profiles show pressure wave propagation and variations in wave height at 312 and 943 seconds, respectively. The axes and colour scale are the same as in Panel A.
As with the previous example, the shallower depth has changed how the pressure above the rupture appears compared with the deeper rupture location § 4.1.2. In fact, this depth has resulted in a dynamic pressure profile that is almost an order of magnitude smaller (note the colour bar in figure 6 compared with figure 2–5). The decreased pressure means the surface-gravity waves dominate the wave profile.
Compared with § 4.2.1, we see that the change in the rupture speed has resulted in the propagation of pressure waves having a shorter wavelength in the
$x$
-direction. We believe that this property leads to pressure waves passing more easily over the trench region, reducing the excitation of standing waves within the trench (figure 6
b). The pressure waves reflected off the ridge structure return to the trench and further excite the standing waves, much like the results of constant rupture velocity. Full animations are provided as Movie 5 (a) and Movie 5 (b).
5. Limitation and future scope
Although the effects of static compression and variable bathymetry are included within a general framework of a finite-depth ocean under linearised water-wave theory, there are several strong assumptions that underpin this work. The first is the assumption of a rigid seafloor. The model could be extended to include the coupling between elastic waves and tsunami waves. Secondly, the generation process can be augmented by including the reflection and transmission of seismic waves radiated by the earthquake at the ocean floor boundary. Further extensions could be made; for example, the assumption of a constant sound profile could be modified to include depth dependence, or nonlinear boundary conditions at the surface could be included. However, we believe that, given the small effect of static compression, the variable sound speed would likely have only a minimal effect, and likewise, the amplitudes of the surface waves are small enough to make nonlinear effects (in the deep ocean) small. We note this model is not suitable to predict near-shore run-up.
6. Conclusion
In this paper, we have created a general formulation for the calculation of fluid pressure and free-surface displacement caused by vertical motion of the ocean floor with an arbitrary piecewise constant bathymetry. We include slight ocean compressibility, including static compression, to facilitate the modelling of acoustic gravity waves. The propagation of AGWs is highly dependent on bathymetry, and therefore the generality of our formulation allows for the straightforward implementation and study of different bathymetric features, such as a continental shelf, underwater trenches and mountain ridges, or any combination of the three. We have shown solutions for bathymetric profiles of a continental shelf and an underwater trench followed by a ridge, with an ocean bed rupture moving vertically at either constant speed or parabolic speed. With just the four cases studied in this paper, we have observed a large variance in pressure wave behaviour in the rupture zone, indicating that the rupture depth and velocity profile must be accurately captured for the appropriate conclusion to be drawn from the simulation results. Over all of our solutions, we have observed that propagating internal pressure waves that appear as vertical bands lead to more reflected waves by any bathymetry rise (e.g. like a continental shelf or ridge). The strong dependence of the pressure profiles on the rupture behaviour and the bathymetry begs the question: How much information can be inferred from measured underwater pressure?
This formulation for an arbitrary depth can be used to expand the results of Gomez & Kadri (Reference Gomez and Kadri2021a , Reference Gomez and Kadri b ), who used a simplified model ignoring the effects of gravity, to a broader scope. The simplicity of the physical set-up poses limitations in capturing the effects of ocean bed properties like elasticity, porosity and viscoelasticity, and the effect of variable sound speed across the water column. However, the computation accuracy of the present method motivates us to incorporate these effects in future research and provide a more holistic approach.
Same as figure 2 at
$t=675.25$
for three difference choices of
$N$
. (a)
$N=15$
, (b)
$N=30$
and (c)
$N=60$
.
Figure 7. Long description
Panel A: A heatmap plot displays pressure wave propagation at 675.25 seconds for N = 15. The x-axis ranges from -150 to 150 kilometres, and the z-axis ranges from 0 to -3.5 kilometres. The colour scale indicates pressure in Pascals, with blue representing lower pressure and yellow representing higher pressure. The surface wave profile above the heatmap plot shows variations in wave height over the same x-axis range. Panel B and C: Equivalent heatmap plot and surface wave profiles show pressure wave propagation and variations in wave height for N = 30 and N = 60, respectively. The axes and colour scale are the same as in Panel A.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2026.11821.
Acknowledgements
S.D. acknowledges the University of Newcastle, Australia, for hosting him as SIRE Fellow, during which the work was initiated.
Funding
This work was partially supported by the Science and Engineering Research Board (SERB), Govt. of India, through the SERB International Research Experience (SIRE) fellowship (Sanction No. SIR/2022/000994) awarded to S.D.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
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.
Appendix. Convergence
To make a remark about the convergence of the mode matching scheme, we present here the solution to the continental shelf bathymetry in figure 7 using three different series truncation values
$N=15$
,
$N=30$
(used for the rest of the paper), and
$N=60$
. We can see that using only fifteen AGW modes does not give a smooth solution near the rupture location; however,
$N=30$
and
$N=60$
give smooth, visually identical results. We note that away from the rupture location, the
$N=15$
solution is very similar to the other solutions.
a0
0
0
−10
−10
t=675.25
N
N=15
N=30
N=60
