On the propagation of acoustic–gravity waves due to a slender rupture in an elastic seabed

Abstract The propagation of waves from a vertical uplift of a slender rectangular fault in a sea of constant depth is discussed, accounting for water compressibility, gravity and seabed elasticity. The compressed water column results in the generation of acoustic–gravity waves that travel at the speed of sound in water. Acoustic–gravity waves are found to terminate after a finite time, with the decay time most influenced by seabed rigidity, which is in contrast to the rigid stationary-phase model where signals persist indefinitely. At certain frequencies acoustic–gravity waves couple with the elastic seabed and travel at the shear velocity (speed of sound in an elastic solid). Improved estimates of the critical frequencies are derived. Moreover, besides the usual tsunami, a second – very small amplitude – surface wave mode travelling at the speed of sound arises under certain frequencies. We derive the cut-off frequency for this mode. The acoustic modes possess a frequency spectrum which depends on the time evolution and spatial properties of the rupture. We find that appropriate filtering of the acoustic–gravity wave signal can reveal characteristic peaks that encode information on the fault's geometry and dynamics.


Introduction
Tsunamis are long-wavelength surface-gravity waves that can cause significant loss of life and damage to property.Recent examples include 2011 Tohoku Oki; 2018 Sulawesi; Palu, Tonga 2022; and the deadliest event so far recorded, the 2004 Sumatra tsunami.Given that tsunamigenic events cannot yet be predicted ahead of time, present-day efforts are directed towards providing early warning systems to mitigate the worst consequences of facing an oncoming tsunami.Existing early warning systems make use of data received from DART (deep ocean assessment and reporting of tsunamis) buoys and seismic measurements.The DART buoys are capable of assessing a tsunami threat, but under standard circumstances, may not leave much time for early warning.Seismic sources can provide information on earthquake size, epicentre etc., but do not comment on possible tsunami characteristics.Tsunamis can be generated by any process that displaces a large volume of water such as volcanic activity, landslides, meteorite strikes and undersea earthquakes.Taking the undersea earthquake as an example, the vertical uplift in the seabed (rupture) elevates and slightly compresses the water column above the rupture zone.Under the restoring force of gravity, fast-moving surface-gravity waves, the tsunami, propagate away from the epicentre at relatively high speed, e.g.200 m s −1 in 4000 m deep water.The compression of the water column above the rupture zone generates low-frequency sound waves, known as acoustic-gravity waves, that propagate at much higher speed ( 1450 m s −1 ) and are able to outrun the tsunami in the far field.Acoustic-gravity waves can couple to the elastic seabed and double their speed (Kadri 2019) and in the solid medium of the seabed compression P waves and shear S waves can propagate at speeds of 6800 m s −1 and 3550 m s −1 , respectively (Dziewonshi & Anderson 1981).However the acoustic-gravity waves in the liquid layer are directly linked to the effective uplift of the rupture zone and thus encode information regarding the rupture geometry and dynamics.In this way the acoustic-gravity waves not only provide a possible means of early detection, but also are a reliable source of information.This information can be extracted, and utilised via an inverse process to determine rupture characteristics (Mei & Kadri 2018;Gomez & Kadri 2021).
Previous work studied the possibility of using acoustic-gravity waves as early warning signals (Yamamoto 1982;Nosov 1999;Stiassnie 2010;Renzi & Dias 2014;Kadri 2015Kadri , 2016)).In these studies the rupture zone has been modelled in a variety of ways (infinite strip, oscillating block, cylinder) and the seabed is normally considered rigid.However, for those cases where the shape of the rupture can be modelled as a rectangular slender body (table 1 of Mei & Kadri 2018) multiple-scales analysis can be applied to take advantage of the differing length scales of the fault.Within Mei & Kadri (2018) the primary focus was on pure acoustic waves, ignoring gravitational effects, whereas a later work (Williams, Kadri & Abdolali 2021) extended the derivations to include gravity and therefore captured the behaviour of the surface-gravity waves in addition to the acoustic-gravity waves.Moreover, the slender-fault model was applied to more complex, multi-fault scenarios, such as the 2016 Kaikoura earthquake (Hamling et al. 2016), via linear superposition -but still with a rigid seabed.However, the elastic properties of the solid medium should not always be ignored (Abdolali, Kadri & Kirby 2019;Kadri 2019) since the water compressibility and seabed elasticity affect the phase speed of surface waves, and thus the arrival times of transoceanic tsunamis (Abdolali et al. 2019).A complementary work by Eyov et al. (2013) investigated the consequences of imposing an elastic seabed as support for a liquid layer residing in a gravitational field upon the form of the dispersion relation.The inclusion of an elastic seabed allows acoustic-gravity waves propagating in shallower water before dissipating into the elastic medium.In stark contrast to the rigid-seabed case the first acoustic mode is able to propagate as a Scholte wave to the shore, where it turns into a Rayleigh wave (Eyov et al. 2013).Note that no rupture was considered in Eyov et al. (2013).The primary objective of this paper is to combine the ground movement of rectangular slender faults with an elastic seabed, in order to study the contribution of elasticity to the propagation of both acoustic-gravity waves and surface waves.The results of this paper should help fill a gap in the literature identified within Renzi (2017) in which the authors claim solutions for acoustic-gravity waves produced by disturbances over an elastic seabed in three dimensions are missing from the literature.
Another difference encountered when studying the elastic case is that rather than the single surface-gravity wave (found in rigid seabed analysis) there is now the possibility of two surface-gravity waves (Eyov 2013).There is the usual tsunami (referred to there as mode 01), and another mode of much smaller amplitude which does not propagate for all frequencies (referred to here as mode 00 (Eyov 2013)).Other secondary objectives include deriving improved estimates for the acoustic-gravity wave critical frequencies, estimating the cut-off frequency for the second surface wave (mode 00), and presentation of a method for rapid calculation of approximate phase velocity curves.We ignore terms of second order and higher (i.e.nonlinear terms) since the free-surface displacements are small in comparison with the water depth (Yoon 2002), and also small in comparison with the wavelengths considered (Michele & Renzi 2020).In addition the gravity term that is present in the full wave equation is omitted because its contribution is small (see figure 2 of Abdolali et al. 2019).
This paper comprises seven main sections.The mathematical formulation combining ground movement with elasticity is found in § 2, with its solution in § 3. Section 4 presents improved approximations for the acoustic-gravity wave cut-off frequencies, and an estimate of the cut-off frequency for the mode 00 surface wave.Section 5 proposes a method for fast calculation of approximate phase velocity curves which does not necessitate solution of the dispersion relation at each data point.Section 6 links the developed theory to numerical results obtained from both synthetic stimulus and real data from hydrophones and DART buoys.The paper concludes with a discussion/summary in § 7.

Formulation
The water layer is considered inviscid, homogeneous, of constant depth h, residing in a gravitational field of constant acceleration g = 9.81 m s −1 .The water layer is assumed unbounded in x and y and is supported by an infinitely deep elastic half-space.The origin of coordinates is taken at the unperturbed free surface directly above the centroid of the slender fault, with the z axis pointing vertically upwards.Assuming irrotational flow, the problem is expressed in terms of a velocity potential function for the liquid φ l , along with a dilation potential φ s and rotation potential Ψ for the solid layer (the subscript s is omitted from Ψ since it only exists in the solid).As in Eyov et al. (2013) we make use of linearised, irrotational flow for the liquid and linear elasticity for the solid.A representation of the flow domain is given in figure 1(a) with a top view of the slender fault in figure 1(b).With i, j, k as unit vectors the velocity in the liquid is given by The solid displacements are then The potentials are governed by three wave equations.In the liquid region: where C l is the speed of sound in water.In the solid region: where C p and C s are the pressure and shear wave velocities respectively: where λ, μ are Lamé constants and ρ s is the density of the solid.At the free surface we have the combined kinematic and dynamic boundary condition: (2.8) In addition, there are four boundary conditions for the seabed.The first of these ensures the vertical component of velocity in the liquid matches that of the solid.The component w s is the vertical component of the seabed motion when there is no rupture (as studied in Eyov et al. (2013)) and as such is small (however, ∂w s /∂t may not be).The magnitude of w s ranges from 10 −6 m for microseisms to 10 −2 m for severe earthquakes (Eyov et al. 2013). (2.9) The definition of W(x, y, t) closely follows that in (3.2a,b) of Mei & Kadri (2018) and describes the motion of the rupture: The duration of the rupture is 2T, the slender fault half-width is b and the slender fault half-length is L. The slenderness parameter is then = b/L 1 (see figure 1b).Note that if there is no rupture, i.e.W(x, y, t) = 0, then the boundary condition (2.9) reduces to that of (8a) of Eyov et al. (2013), ẇl = ∂w s /∂t.On the other hand, when the seabed is rigid, w s = 0, and we recover the bottom boundary condition (2.3) of Mei & Kadri (2018), ẇl = W(x, y, t) = −∂φ l /∂z, but this time with a minus sign due to this paper following the sign choices in (1a) and (1b) of Eyov et al. (2013).The next boundary condition states that the axial stress σ zz is equal in magnitude but of opposite direction to the liquid pressure at the seabed: (2.12) The remaining two boundary conditions define no shear on the seabed: The dynamic pressure and surface elevation are obtained from We also require φ l , φ s , Ψ and all derivatives to decay to zero as x, y, t → ±∞, z → −∞.

Solution
We introduce multiple-scale coordinates following Mei & Kadri (2018): The wave equations (2.4)-(2.6)can then be written as X, Y, z, t), with similar expressions for φ s and Ψ .
Then the perturbation equations at O( 0 ) describe the two-dimensional problem of an infinitely long slender fault: The fault motion, elastic properties and elastic dispersion relation are all captured at O( 0 ).Thus, the O( 2 ) boundary conditions for the liquid layer are those for rigid seabed and no fault motion: (3.12)
To arrive at a trial solution for (3.14) we choose φs0 (z) = D 1 e qz , because φs0 (z) → 0 as z → −∞ implies no terms involving e −qz .By a similar argument Ψ 0 (z) = D 2 e sz j.
Applying the double Fourier transforms to the boundary condition at z = 0 (leading-order term): Again, after application of both Fourier transforms we arrive at (3.28)

Transformed governing equations
We assemble terms and drop the zero subscript for ease of notation (but remembering these are leading-order terms).Also, we drop the double overbar, again remembering that these are the terms after the double Fourier transforms: (3.30a-c) At z = 0 we have the (transformed) boundary condition for the liquid surface: Then, at z = −h we have four (transformed) boundary conditions for the seabed: 3.3.Form for potentials Taking Φ l (z) = E 1 cos(irz) + E 2 sin(irz), we substitute into the boundary condition at z = 0 to arrive at (3.36) in agreement with (14a) of Eyov et al. (2013).Also, we take Φ s (z) = D 1 e qz + D1 e −qz and Ψ (z) = ψ y j with ψ y = D 2 e sz + D2 e −sz , but note that, in order to obtain physical solutions in which solid displacements decrease with depth, we must have D1 = D2 = 0 and s, q ∈ R ≥0 .If this were not the case then displacements would oscillate or increase with depth (Eyov 2013).
Applying boundary condition σ xz = 0 (3.34) at z = −h: (3.37) Applying boundary condition (3.32) at z = −h: Applying boundary condition (3.33) at z = −h: Since (3.38) and (3.39) are essentially a pair of simultaneous equations in unknowns E 1 and D 1 they can be solved in this case, resulting in with ) and Setting H 2 = 0 and rearranging yields which is the dispersion relation (17) of Eyov et al. (2013).The zeros of H 2 (i.e.dispersion relation solutions) locate the poles for the residue calculations that come later.Therefore, we have Setting q = s = 0 reduces to the rigid case where (3.1) of which is in agreement with (3.1) of Williams et al. (2021) (note that the sign difference is due to the definition of the velocity potential).

Inverse Fourier transforms
The leading-order potentials are retrieved by applying the inverse Fourier transforms as where I 1 , I 2 , I 3 are the k integrals: (3.51a,b) In each case the integrand has poles at the zeros of H 2 , i.e. whenever the dispersion relation (3.44) is satisfied.We substitute out r, q and s to make H 2 purely a function of k.Then values for I 1 , I 2 , I 3 can be calculated from the residues.
Figure 2 shows the various zones where r, q and s take on real and imaginary values.There are zones corresponding to surface waves and acoustic-gravity waves.
Zones possible according to r, q, s being real or imaginary for the case ω = 2π, C l = 1450 m s −1 , C s = 3550 m s −1 , C p = 6300 m s −1 .Zone 1 (orange) has r, q, s ∈ R and corresponds to surface-gravity waves.Zone 2 (green) has r ∈ iR, with q, s ∈ R and corresponds to acoustic-gravity waves.
The remaining zones near k = 0 (grey) are not physical solutions.The points where r, s, q transition real imaginary are designated ±k r = ±0.00433(black circles), ±k s = ±0.00177(red circles) and ±k q = ±0.00099(blue circles) respectively.
The remaining zones close to k = 0 are not physical solutions, since imaginary values taken on by q and/or s would imply oscillations at infinite depth in the elastic medium.Moreover, q and s have to be real and non-negative, otherwise oscillations would increase with increasing depth into the elastic medium.Examination of I 1 , I 2 , I 3 indicates that possible poles might also exist at k = 0 and when k 2 + s 2 = 0.When k = 0 the sin(kb) term in the numerator (from H 1 and H 3 ) ensures a factor of b is reached in the limit k → 0, so k = 0 is a removable singularity.For the case k 2 + s 2 = 0 there is a possible pole when k = ω/ √ 2C s , but this pole lies in the unphysical zone of figure 2. From Eyov et al. (2013) we have that s = 0 (at k s ) represents a point where the energy spreads out over the whole solid depth.At that point the wave amplitude vanishes and so propagation ceases.
When r ⇒ r 0m with m = 0, 1, which corresponds to surface waves.There are two possible modes for surface waves.
Mode 00 can propagate if ω > ω 00 -the cut-off frequency for this mode.Mode 01 is the usual tsunami.If instead r ⇒ ir n , then acoustic-gravity waves are possible, and up to a maximum value of n = N, after which the evanescent waves exist with wavenumber Λ n : Solutions to the dispersion relation involving acoustic-gravity waves for the case ω = 2π occur between k s = 0.00177 and k r = 0.00433.They are marked with blue diamonds in figure 3. Considering the liquid terms first, we break φ l0 into the different regions according to varying ω.For r ∈ iR: whereas for r ∈ R: Application of the Rayleigh damping method and contour integration using the residue theorem around a positively oriented simple closed curve as per Mei & Kadri (2018) results in with The derivative terms are given in Appendix A.
In support of the validity of the integration process figure 4 shows a plot of 1/|H 2 | in the complex plane when H 2 = H 2 (k) and k is allowed to take on complex values.Cross-sections through the real and imaginary axes appear in figures 5(a) and 5(b) respectively.The poles of the function lie on the real axis, whereas the zeros lie on the imaginary axis.If the range of the plots were to be extended then the function decays to zero everywhere.As empirical evidence for the validity of the integration, when the calculations are complete, we find good agreement with existing synthetic and real data plots for both acoustic-gravity waves and surface waves (e.g.see figures 13 and 17).
In the case where r ⇒ r 0m , k ⇒ k 0m with r 0m a real number and m = 0, 1 then there may exist two possibilities for surface waves: gr 00 sin(irz) e ik 00 x , (3.62) with Acoustic-gravity waves from slender rupture in elastic seabed The derivative term is again to be found in Appendix A. Using the substitutions the dispersion relation (3.44) can be written in terms of r and ω alone, and in this case the condition for the existence of the 00th mode for a particular ω is d dr left-hand side of (3.44) > d dr right-hand side of (3.44) . (3.66) The expressions for the velocity potential in the liquid layer can be further reduced to (3.67) Returning to the expression for the displacement potential in the solid given by (3.68) and following a similar procedure to that of the liquid velocity potential case, we arrive at with The derivative terms evaluated at k = k n , k = iΛ n and k = k 0m remain as before.In a similar way the rotation potential can be written as which becomes k=k 00 e h(s 00 −q 00 )+s 00 z e ik 00 x (3.72) 3.5.Long-range modulation: liquid layer We introduce unknown envelope factors for the liquid layer A ± n (X, Y) and A ± 0m (X, Y): (3.73) The initial conditions are given by The waves are required to vanish far away from, and be symmetric about, the central axis (Mei & Kadri 2018): Proceeding with acoustic modes, taking the time Fourier transform of (3.8) and separating φl2 into three ranges yields (3.76) In the range From this point the solution proceeds in an analogous way to that derived by Mei & Kadri (2018): gr sin(irz) e ik n x . (3.78) Assuming φ+ l2 has solutions in the form N n=1 ξ + n (ω, z) e ik n x , then substituting into (3.78)gives Equating coefficients gives where where gr n sin(r n z). (3.83) The ground motion is captured at O( 0 ), so the boundary conditions on φ+ l2 (and therefore Here F n (z) is a solution to the boundary value problem: ) A similar process could be carried out for surface waves.The next step is to extract the Schrödinger equation from (3.82) by applying the following Green's identity over the range −h ≤ z ≤ 0: As in Mei & Kadri (2018), the result is the Schrödinger equation for the two-dimensional evolution of the envelope factors: Acoustic-gravity waves from slender rupture in elastic seabed here as where C(z) and S(z) are Fresnel integrals.A similar process beginning at (3.76) can be applied to derive the expressions for the surface wave mode 01 and mode 00.Finally, the pressure obtained from the velocity potential (3.67) in the liquid (propagating parts) inclusive of envelope factors is given by Similarly the surface elevation is given by (3.93)

Improved critical frequency approximations
In a practical application of (3.92) and (3.93) numerical solutions approximate the integrals over a finite range, and so knowledge of the critical frequencies ω n and ω 00 is essential.The critical frequencies ω n represent the cut-off for acoustic-gravity wave mode numbers n ≥ 2, and ω 00 is the cut-off for the surface wave mode 00.The first acoustic-gravity wave mode does not have a cut-off frequency (see figure 6b).An approximation for ω n exists in the form of (4.1) (Eyov et al. 2013), but this approximation is based upon the location of the vertical asymptotes found in the dispersion relation plots, an example of which is shown in figure 6(a).This approximation -although compact and easy to use -is not as accurate as it might be.The following subsections construct a more accurate approximation for ω n (albeit more complicated) and an approximation to the cut-off frequency of surface wave mode 00 based on the gradient condition (3.66).progressive mode (n = 1) for an elastic seabed is a Scholte wave, which propagates all the way to the shore, where it turns into a Rayleigh wave.From (30) of Eyov et al. (2013), the critical frequency for a particular depth is given by

Acoustic-gravity waves
into (4.2) to obtain a function of r alone, followed by another substitution r ⇒ ir to retrieve the acoustic-gravity wave solutions. Let , which is the s = 0 condition for the termination of progressive modes, and then let r = (n − 3 2 )π + δ(n), where δ(n) represents a small (mode-dependent) positive offset away from the vertical asymptotes located at r = (n − 3 2 )π.So the desired δ(n) is the r separation between the blue and red circles in figure 6(a).Ignoring terms of O(δ(n) 2 ) an approximation of the dispersion relation (see Appendix B) can be written as where the coefficients a n , b n , c n , d n are given in Appendix C. Then using the approximation − cot(δ(n)) −1/δ(n) for small δ(n), (4.5) can be put into quadratic form: This can be solved for δ(n), and then the value of rn and the critical frequency ω n can be obtained from To determine how well δ(n) predicts the offset a comparison was made between the approximate value of ω n calculated using (4.7) and that found by using the dispersion relation The comparison was carried out for depths ranging from 500 to 8000 m and all available modes.The maximum error occurred in the second mode at a depth of 8000 m, but was still less than 0.1 % (see figure 7).In § 5 the results for rn and δ(n) are used to construct approximate phase velocity curves.

Surface wave
The surface gravity mode 00 does not exist for all frequencies in the elastic seabed case and never exists for a rigid seabed.The cut-off condition for mode 00 is given by (3.66) which can be solved numerically for a given depth.Alternatively a good approximation can be obtained by seeking the frequency at which the gradient of the left-hand side of the dimensionless dispersion relation in (4.2) is equal to the gradient of the right-hand side.This occurs for small (ultimately zero) r.In this case we make the approximation tanh(r) r.Now we differentiate (4.2) with respect to r.Then we express the result as a series in r2 to arrive at (4.9) In the limit r → 0, (4.9) can be written as the quadratic ω2 − A ω − 1 = 0, with (4.10) Taking the positive root of the quadratic gives the approximate cut-off frequency which we name Ω00 .The dimensional form can be recovered from Ω 00 = Ω00 g/h.A workable approximation to the cut-off frequency can be obtained by taking A. We call this approximation A 00 .2a) or (ii) fix a constant depth and plot phase velocity versus frequency (as in this paper).In either case for each data point on every curve the dispersion relation has to be solved numerically which can be time-consuming.Also, care has to be taken to ensure solutions are valid.Here we present an alternative method for quickly plotting an approximate version of the elastic seabed phase velocity curves.In the following, variables with a tilde are made dimensionless according to (4.3).The method is based around the tanh −1 function which is manipulated in the following ways.
• Scale along the horizontal r axis (the independent variable) so as to fit the range (n − 3 2 )π . . .(n − 1 2 )π, with n being the mode number: (5.1) • Shift up the vertical axis so that at centre range r = (n − 1)π the value is αCs, i.e. the Rayleigh wave phase velocity where the first acoustic mode intersects the vertical axis (see figure 6b).The value of α = 0.922231 is taken from (5.22) of Eyov (2013): (5.2) • Next stretch the plot along the vertical axis by a factor κ(n) so that the curve meets the shear velocity C s at the critical value rn determined from (4.7): where and δ(n) is the critical offset calculated via the procedure described in § 4. • Include a multiplicative factor Ỹ(r, n) to ensure that the curve has its region of rapid descent shifted away from the (n − 1 2 )π asymptote to better align with the 'reference' phase velocity curves derived from the dispersion relation, and so help minimise errors.The function ṽ(r, n) so obtained is the generating function from which all the phase velocity curves are derived: where (5.6) The resulting curve for the case n = 1 is shown in figure 8. Curves for higher modes are obtained by shifting the n = 1 case by appropriate multiples of π along the positive r axis.The variable r ranges over the interval rn ≤ r ≤ (r * − ε) with 0 < ε 1, defined by 2 )π and r * is such that ṽ(r * , n) = Cl .This represents the phase velocity asymptotically approaching C l with increasing frequency (all modes).• Take the generating function for each mode and translate, so that the known point ( ωn , Cs ) → (0, 0).This is the black curve t(r, n) in figure 9: (5.7) where zn is a fixed complex number representing the known cut-off point ( ωn , Cs ).• Apply the shearing function S(ã, n) to distort the black curve into the appropriate shape for the mode considered -the coloured curves zt n in figure 9: (5.9) A plot of S(ã, n) is shown in figure 10(b).The function w(ã, n) is derived from the phase velocity curves for the rigid seabed case (figure 10a) by inverting (5.10) to give ω in terms of rigid seabed phase velocity ṽr : (5.10) The expressions for kn and ωrn appearing in (5.10) are from (3.9) and (3.10) of Mei & Kadri (2018) here made dimensionless.After performing the inversion (5.12) The function ã is a measure of the vertical drop from the constant line Cs down to the phase velocity curves ṽr (see figure 10a).In order to apply S(ã, n) to the generating function define ã = C s − ṽ(r, n) so now ã represents the vertical drop from the constant line Cs down to the phase velocity curves ṽ(r, n) (see figure 8).• Add zn to translate back: (5.13) • Finally re-scale to obtain the desired phase velocity curves -solid black trace in figure 11: (5.14) The solid black trace of figure 11 is a complex plot with real part representing the angular frequency ω and imaginary part representing the elastic seabed phase velocity v e .The dashed curves of figure 11 are those obtained by numerically solving the dispersion relation (3.44).To quantify the errors between the phase velocity obtained using the shearing method and that obtained by solving the dispersion relation, use (5.15) The maximum error occurs for the first mode and errors decrease with increasing frequency and increasing mode number (figure 12).There is some freedom in the expression for Ỹ(r, n) which could potentially reduce the errors a little by carefully replacing the π/18 term with an alternative value derived from some error minimisation technique (e.g.minimax approximation (Powell 1996)), which we did not pursue further here.

Numerical results
6.1.Acoustic-gravity waves Figure 13 compares the first acoustic mode for the elastic case (3.92) and rigid case ((3.22) of Williams et al. (2021)).As in Eyov et al. (2013) the values used for ρ l , ρ s , C l , C s and C p are average values taken from Dziewonshi & Anderson (1981).One stark difference between the two cases is that the signal terminates after some time in the elastic case, whereas it continues indefinitely in the rigid case.Another difference is the presence of signal at times earlier than the main pulse in the elastic case, but no signal at all in the rigid stationary phase model.The phase velocity curves for the elastic case (figure 11) indicate that frequencies close to the critical frequency for each mode receive a boost in phase velocity enabling signals to propagate faster.For these frequencies speeds close to C s are achievable.The rigid seabed stationary phase model produces complex numbers for times earlier than x/C l due to a singularity induced by the stationary phase method (Stiassnie 2010).The pressure amplitudes are similar.
A sensitivity analysis was carried out looking at the effects of six parameters on the signal duration (see figure 14).Each parameter was varied individually away from its reference value (table 2) while holding all other parameters at reference.Then the percentage change in pulse duration was divided by the percentage change in the parameter to arrive at the sensitivity value.It was found that the rigidity of the seabed most affected the signal duration.Increasing the Lamé parameters increases C s and C p in accordance with (2.7a,b); the ratio C p /C s was kept constant.Another difference between the elastic seabed case and the rigid seabed case appears when a fast Fourier transform of the signals is examined.The elastic case shows a slight upward shift in the frequency peak (see figure 15).This is in contrast to the slight downward shift found when a viscous compressible sediment layer is overlying the seabed as in Abdolali, Kirby & Bellotti (2015).Whilst investigating the synthetic acoustic-gravity waves, we found that band-pass filtering applied to the signal generated by combining the first 10 modes (figure 16a) revealed some interesting peaks located close to the expected arrival time for a phase velocity of C l (figure 16b).There are four peaks of particular interest labelled 1, 2, 3 and 4. The presence of the peaks is a consequence of the fault's geometry and motion, and is not related to the rigidity of the seabed, since the peaks are also present under rigid seabed conditions, and even when the signal considered was purely acoustic, as in Mei & Kadri (2018).The time spacing between pairs of peaks responds to changes in either fault half-width b or rupture duration τ = 2T in a linear fashion, so that details of the fault's geometry and dynamics are encoded in the acoustic-gravity waves.Time t 1 between peak numbers 1 and 2 (or 3 and 4) is exactly the rupture duration, and t 2 between peaks 1 and 3 (or 2 and 4) is linearly related to the fault half-width through t 2 = 2b/C l .When the slender fault begins to move, peaks 1 and 3 are generated at the edges of the slender fault and begin to propagate.The time separation between these peaks is explained as the time required for a wave travelling at speed C l to cross a fault width of 2b.At the end of the fault's motion, after τ seconds, the second pair of peaks (2 and 4) is generated and propagates away -also separated in time by t 2 .The resultant waveform as would be seen in the far field is a collection of four peaks.The amplitude of the peaks depends linearly on the uplift velocity W 0 .The timings between peaks agree well with the values given in table 2. Let subscripts 1, 2, 3 and 4 represent the peaks denoted by the numbers 1, 2, 3 and 4 respectively.Then, t 12 = 10.97 s, t 34 = 8.9 s, t 13 = 57.67 s and t 24 = 55.46 s.The times t 12 and t 34 represent τ , the uplift time, which (from table 2) is actually 10 s.The times t 13 and t 24 are the transit times for an acoustic signal to cross a fault width 2b, which, again from table 2, is actually 55.17 s.The information that could be extracted from the timings embedded in the acoustic modes would be of interest to the inverse process that reconstructs fault parameters from received signals (Gomez & Kadri 2021).In § 6.3 an actual hydrophone recording made during the Samoa 2009 event is filtered to reveal the four peaks encoded within it.Nosov (1999)   of the bottom displacement -which is referred to as tsunami's voice.This is exactly what we find encoded into the characteristic (four) peaks.
6.2.Surface waves Consider now the surface waves generated by the single slender fault with parameters as per table 2. The equations generating the surface waves are (3.93) for the elastic case and (3.23) of Williams et al. (2021) for the rigid case.At a depth of h = 4000 m there is little difference between elastic and rigid cases (figure 17a), but at a depth of h = 1000 m differences are more apparent (figure 17b).In deeper water the surface wave is almost unaffected by the elasticity of the seabed (Eyov 2013).This surface wave is the main tsunami, i.e. mode 01.
When the seabed is elastic the possibility of a second surface wave arises.This wave does not exist for all frequencies and never exists in the rigid case (Eyov 2013).The gradient condition (3.66) has to be satisfied before mode 00 can propagate (see figure 18).The mode 00 surface wave has a phase velocity C l , and a negligible amplitude, of the order of micrometres.A plot of mode 00 under the conditions of table 2 can be found in figure 19.In the plot there are four distinct peaks numbered 1 to 4. These peaks in the mode 00 surface wave correspond to the peaks numbered 1, 2, 3 and 4 found in the acoustic signal.The assumption of a rectangular fault moving at a uniform speed results in symmetric peaks, whereas in reality the motion is much more complicated thus upsetting the symmetry of the peaks seen in figure 19.

Hydrophone recordings
The theory developed in this paper leading to the equations for pressure (3.92) and surface elevation (3.93) is linear.Therefore, as in Williams et al. (2021), more complicated multi-fault scenarios can be constructed from single slender fault solutions by linear superposition, given that the parameters for each individual fault are known.To contrast the case of an elastic seabed with that of a rigid seabed we revisit the Sumatra 2004 earthquake discussed in § 5.1.2 of Williams et al. (2021).The geographical area considered here ranges over [70 • to 100 • ] E longitude and [−15 • to 20 • ] N latitude.This curved patch on the (idealised) spherical Earth is mapped to flat x, y coordinates.The conversion factor of metres per degree is fixed for the latitude (y) direction, but the metres per degree in the longitudinal (x) direction varies with latitude, being maximum at the equator and decreasing as the poles are approached.To simplify the calculations an average value for metres per degree longitude was used and the area considered was kept reasonably small.The fault centroids are marked by black stars in figures 20(a) and 20(b) and all faults are contained within the masked off 'earthquake zone'.The purpose of the earthquake zone was to avoid pressure calculations too close to the faults.The location of hydrophone H08N is marked with a red star.Figure 20  the earthquake zone.The elastic seabed could therefore be considered more physically realistic.
Figure 21 compares the predictions made by the elastic seabed model against recorded data for the Sumatra event derived from the southern (H08S1) hydrophone and the the acoustic-gravity waves in the elastic model, the model is unable to predict the P-wave portion of the hydrophone signal.Modifications to the existing model or maybe a new model would be needed to capture this behaviour.Between the green and blue lines the elastic model predicts acoustic-gravity waves that can travel at phase speeds close to C s for some frequencies.The hydrophones show weak signal in this region, possibly due to the filtering effect of the hydrophone's response.After the blue line the elastic model predicts a signal that is close in amplitude to that of the hydrophone recordings, but decays more slowly -possibly due to a lack of dissipation included in the model.However, the signal duration is at least finite in the elastic case.There are processes missing from the elastic model (varying bathymetry, reflection, refraction, dissipation etc.) so an exact match is not expected.Examination of figure 21 shows the arrival times for P-waves, S-waves and the main acoustic-gravity wave pulse (travelling at a phase speed of C l ) are consistent with our assumptions of constant water density, constant speed of sound in liquid and constant speed of propagation in solid.The leading pulse seen in the hydrophone recordings is primarily made up of lower-frequency components.To show this, a band-pass filter was applied to the hydrophone recording at H08S.The filter eliminates most of the frequency components below 3 Hz and has largely flattened the leading pulse of the H08S signal (figure 24).In order to enhance the detection of acoustic signal between the red and blue lines, ultralow-frequency hydrophones should be used.
To demonstrate the extraction of fault timing and geometry from acoustic-gravity wave signals a band-pass filter was applied to the data obtained from hydrophone H11 located at Wake Island during the Samoa 2009 event (data supplied by CTBTO).The timings for the peaks are t 12 = 15.46,t 34 = 27.75, t 13 = 29.89 and t 24 = 42.18s (figure 25).Note that the time axis in figure 25 does not begin from the start of the rupture.The time axis in this case represents an 1800 s window around the main hydrophone signal.This is not a concern here since only time differences ( t values) are required.Unlike in the synthetic case this event is asymmetric in the sense that the timings suggest either a trapezoidal rupture geometry or non-uniform uplift velocity (or both).The tsunami profile is affected by seabed elasticity in shallower water.The inclusion of elasticity induces a decay into the acoustic-gravity wave signals so that the signals terminate after some finite time, unlike the rigid, stationary phase model.From the parameters studied, we find that the signal duration is most affected by the seabed rigidity, with duration increasing alongside rigidity until the totally rigid condition is achieved, at which point the signals persist indefinitely.Thus the inclusion of elasticity helps facilitate a more realistic representation of the pressure field.When the seabed is elastic there exists the possibility of two surface waves.The first (mode 01) is the usual tsunami, but the second (mode 00) is an interesting mode which does not propagate for all frequencies in the elastic case, and never exists in the rigid case.Linear relationships between mode 00 timing of signal peaks and the fault parameters b (half-width) and τ (rupture duration) are found and explained.There is also a linear relationship between uplift velocity and mode 00 amplitude.Information on the fault geometry and timing is encoded into the mode 00 surface wave, and is also found to be imprinted into the acoustic-gravity wave signals as well.With appropriate filtering it is possible to extract this information from the acoustic-gravity wave signal (at least in some instances), which would be helpful in solving the inverse problem of deriving fault properties from acoustic/seismic information.Additionally, an improved estimate of the critical cut-off frequency for acoustic modes n ≥ 2 is presented, which is then used in a new method for calculating approximate phase velocity curves which does not rely on solving the dispersion relation (3.44) at each point.The facility to quickly produce approximate phase velocity curves may help in reducing the computational burden in real-time analysis.In a side-by-side comparison of the shearing method against the dispersion solving method, the shearing method was found to be approximately twice as fast.The comparison was run on the same computer (Intel(R) i9 CPU, 3.60 GHz, 128 GB RAM) and used the same software (Maple) to produce phase velocity curves for 16 modes with ω = 20 m s −1 and depth h = 4000 m.An approximation for the mode 00 surface wave cut-off frequency is also derived.As in many previous studies, the model developed here has a constant water depth assumption, so while the model can determine the tsunami properties for deep water, it may fail for varying bathymetry.changes in bathymetry without computation of the entire three-dimensional domain.For slowly varying bathymetry (i.e.mild slopes where ∇h(x, y, t) kh) there already exist techniques in the form of the depth-integrated equations (Sammarco et al. 2013;Abdolali et al. 2015;Renzi 2017).In the conclusion to Renzi (2017) the authors remark that models of tsunamigenic events over an elastic seabed do not appear in the literature to date.This topic is addressed and solved (at least for constant depth) in § 3 of this paper.
In the study of the bottom pressure field for the Sumatra 2004 event covered in § 6, a curved patch of the Earth's surface was mapped to a flat x, y plane.An interesting extension to this work could be to move the perspective of the study into a more global viewpoint by use of spherical coordinates.In that way far-field predictions may become more accurate.

956Figure 1 .
Figure 1.Representation of the flow domain.(a) Cross-section through the x, z plane.Water depth is h, surface elevation is η(x, y, t), liquid velocity potential is φ l , solid dilation potential is φ s and solid rotation potential is Ψ .Densities in the liquid and solid medium are ρ l and ρ s respectively.(b) Top view of slender fault.

Figure 3 .
Figure 3. Acoustic-gravity wave solutions to the dispersion relation are located at the intersections of dashed and solid curves (blue diamonds) for ω = 2π and depth h = 4000 m.Dashed curve is left-hand side (LHS) of (3.44) and solid curve is right-hand side (RHS) of (3.44) when r ∈ iR.

Figure 4 .Figure 5
Figure 4. Plot of 1/|H 2 | in the complex plane when H 2 = H 2 (k) and k is allowed to take on complex values.The angular frequency in this case is ω = 2π as in figure3.

Figure 6 .
Figure 6.Approximate critical values from (4.1) (red circles) and actual critical values (blue circles).(a) Dispersion relation plot for h = 4000 m.Red circle marks vertical asymptote.Blue circle marks r 2 , the actual cut-off for mode 2. Dashed trace, left-hand side (LHS) of equation (3.44); solid trace, right-hand side (RHS) of equation (3.44).(b) Phase velocity curves for first four modes at constant depth of h = 4000 m.Dotted line is C s = 3550 m s −1 .
.1) which is a good approximation to the actual critical frequency, though it is based on the location of the vertical asymptotes in the dispersion relation plot -location of red circle in figure6(a).For accuracy, we require a better approximation to the actual intersection of the two curves in the dispersion relation plot -blue circle in figure6(a).We begin with 956 A6-20 https://doi.org/10.1017/jfm.2022.

956Figure 7 .
Figure 7. Percentage error for approximate critical frequencies ω n from (4.7).Depths range between 500 m (lower error bound) and 8000 m (upper error bound) -all available modes.
.11) 956 A6-24 https://doi.org/10.1017/jfm.2022.functionṽ(r, n) for first acoustic-gravity mode (n = 1) with depth h = 2000 m.Other modes are derived by shifting the horizontal axis through (n − 1)π and using the appropriate values for rn and r * .make the change of variable ṽr = Cs − ã, and rename ω to w to arrive at
Figure 14.Response of signal duration when changing parameters.

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

956Figure 18 Figure 19 .
Figure18.Left-hand side (LHS) of dispersion relation (3.44) (dashed trace) and right-hand side (RHS) of dispersion relation (solid trace) when r ∈ R. The frequency is at the point where mode 00 becomes active, ω 6.95 rad s −1 , h = 4000 m (see table1).The top logarithmic plot indicates overall behaviour.The middle and bottom plots provide expanded views.The mode 00 solution in the middle plot is where the solid curve touches the dashed curve 0 ≤ r ≤ 0.1, and the mode 01 solution (the usual tsunami) in the bottom plot is where the solid curve again makes contact with the dashed curve in the descending phase 2.5 ≤ r ≤ 3.

Figure 20 .Figure 20
Figure 20.For caption see next page.

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

Table 1 .
Table1gives values for ω 00 found using a numeric solver and compares with the approximations Ω 00 and A 00 indicating errors for various depths.Comparison of cut-off frequencies obtained from numeric solver (ω 00 ) with approximations from quadratic solution (Ω 00 ) and coarse approximation (A 00 ) for various depths h.
Eyov et al. (2013).1017/jfm.2022.10915.Approximate phase velocity curves: shearing methodWhen plotting phase velocity curves it is typical to choose one or other of the following scenarios: either (i) fix constant frequency ω and plot phase velocity versus depth h as inEyov et al. (2013)(figure

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