Acoustic–gravity waves from multi-fault rupture

Abstract The propagation of wave disturbances from a complex multi-fault submarine earthquake of slender rectangular segments in a sea of constant depth is discussed, accounting for both water compressibility and gravity effects. It is found that including gravity effects the modal envelopes of the modified two-dimensional acoustic waves and the tsunami are governed by the Schrödinger equation. An explicit solution is derived using a multi-fault approach that allows capturing the main peak of the tsunami. Moreover, a linear superposition of the solution allows solving complicated multi-fault ruptures, in particular in the absence of dissipation due to large variations in depth. Consequently, the modulations of acoustic waves due to gravity, and of tsunami due to compressibility, are governed simultaneously and accurately, which is essential for practical applications such as tsunami early warning systems. The results are validated numerically against the mild-slope equation for weakly compressible fluids.


Introduction
Tsunamis have a long history of devastation costing the lives of thousands of people, causing damage to property and posing a continuing threat to thousands of kilometres of shoreline.In the last two decades we have witnessed the deadliest 2004 Sumatra earthquake and tsunami, followed by other major [tsunamigenic] earthquakes, such as the 2011 Tohoku Oki, and more recently, the 2018 Sulawesi and Palu tsunami.These have repeatedly proven that a primary reason for the scale of devastation is the failure to provide a reliable early warning.Current warning systems rely heavily on DART buoys (Deep ocean Assessment & Reporting of Tsunamis) and seismic measurements.Accurate tsunami evaluation from DART buoys may be possible, though depending on particular circumstances there may not be sufficient time for early warning.On the other hand, seismic data provide valuable information on the tectonic movements, earthquake size, and possible epicentre, though with current technology and analysis they fail to assess the tsunami threat.A complimentary approach has been suggested by a number of authors who considered the slight compressibility of the water in their analysis (Miyoshi 1954;Sells 1965;Yamamoto 1982;Nosov 1999;Stiassnie 2010;?;Kadri & Stiassnie 2013b;Cecioni et al. 2014;?;Abdolali et al. 2018).In this approach, attention is focused on acoustic-gravity waves that radiate from submarine earthquakes alongside the tsunami, and propagate through the liquid or elastic layers (Eyov et al. 2013;Kadri 2016).Acoustic-gravity waves are compression waves that reside within the entire water column and can couple with the elastic sea bottom.They carry information about the source at relatively high speeds ranging from the speed of sound in water (1500 ms −1 ), to Rayleigh waves speed in the solid (3200 ms −1 ) that far exceeds the phase speed of the tsunami (200 ms −1 at 4 km water depth), see Eyov et al. (2013).In the solid layer, compression P (pressure) waves and S (shear) waves propagate at about 6,800 ms −1 and 3,900 ms −1 , respectively (e.g.see preliminary earth reference model Table 1 of Dziewonshi & Anderson (1981)).A critical difference between analysing acoustic-gravity waves and P and S waves is that the former, being a compression wave in the liquid layer, is directly associated with the effective vertical uplift.Hence, acoustic-gravity waves can not only act as excellent precursors, but they could also provide vital information on the geometry and dynamics of the effective uplift, which eventually shapes the main characteristics of the tsunami.Note that sea-bottom elasticity can significantly affect the phase speed of acoustic-gravity waves but only in shallow water (Eyov et al. 2013;Kadri 2019).On the other hand, water compressibility and sea-bottom elasticity affect the phase speed of surface gravity waves and should be considered when accurate transoceanic tsunami modelling is sought (Abdolali & Kirby 2017;Abdolali et al. 2019).The peak frequency shift and attenuation of acoustic-gravity waves takes place due to interaction with a visco-elastic sedimentary layer at the sea bottom (Abdolali et al. 2015b;Prestininzi et al. 2016).
The possibility of using the acoustic fore-runner as an early warning signal has long been established (Nosov 1999;Yamamoto 1982;Stiassnie 2010;Kadri & Stiassnie 2013a; ?; Abdolali et al. 2015a;Kadri 2016;Kadri & Akylas 2016;Renzi & Dias 2014).Finite fault models have also been investigated, providing a three dimensional theory of acoustic-gravity waves based on the classical method of the Green's function (Hendin & Stiassnie 2013).However, their utility in providing predictions for acoustic and surface wave behaviour in real time is limited.The limitation arises due to the solution being in integral form which requires partitioning of any shape considered into many small elements, calculating the contribution from each element, then performing a summation to arrive at the total contribution.In the absence of an explicit analytical solution, this proves to be computationally expensive (Mei & Kadri 2018) and, of course, the processing burden escalates with the addition of more complex multi-fault ruptures, as observed in nature (Hamling et al. 2017).An alternative approach was proposed by Mei & Kadri (2018) who considered a slender fault and invoked multiple scales analysis to obtain a closed form analytical solution for the propagating acoustic modes.Improvements in long-range modulation are provided by the introduction of envelope factors involving Fresnel integrals.Moreover, an inverse approach was employed and relations for fault location and rupture parameters were derived.
There are to main objectives of the present work.The first, is to extend the results of Mei & Kadri (2018) to include gravitational effects and multi-fault ruptures.The inclusion of gravitation involves a modification to the surface boundary condition.This modification gives rise to expressions for the gravity wave contribution to bottom pressure, along with the expected acoustic-gravity wave contributions.Evanescent modes are also derived, but later ignored, since their effects in the far field are negligible.Expressions for surface elevation are obtained -broken down into contributions from the surface wave and the acoustic-gravity waves.The form of the governing equations for the envelope factors involved in the long range modulations are found to be identical for both acoustic-gravity waves and the surface wave, i.e. they both obey the Schrödinger equation.The addition of gravitation to the current model may have a beneficial effect on the accuracy achievable in the inverse theory calculations originally discussed in Mei & Kadri (2018), which is an ongoing work.The second objective is to tackle a long standing limitation which arises when applying a stationary phase approximation.The derived explicit solution for the gravity mode (tsunami) is singular at the arrival time which results in overlooking the main peak of the tsunami.To overcome this difficulty, we employ a multi-fault approach, where the original fault is split into stripes.Since each stripe has a different spatio-temporal singularity, the main tsunami amplitude can be reconstructed by superposition principle.
Extension to multi-fault ruptures arises naturally from the linear theory by application of the superposition principle enabling fault systems such as that discussed in Hamling et al. (2017) to be investigated.Two instances of multi-faults are considered here, one based upon the 2011 Tohoku event (detailed data can be found in Abdolali et al. (2017)) and the second is based on the Sumatra 2004 event.
In this paper we explicitly ignore contributions from terms of second order and higher, (i.e.non-linear terms) on wave amplitude over the spatial range of interest.In addition we assume that the evolution of each acoustic mode does not involve mode coupling (Michele & Renzi 2020).This is well justified for acoustic-gravity waves that are the main focus of this work.

Governing equations
The water layer is assumed to be inviscid, homogeneous and of constant depth h.The origin of the Cartesian co-ordinates is located at the sea-bottom, at the centroid of the fault, with the vertical axis z directed vertically upward.Based on irrotational flow, the problem is formulated in terms of the velocity potential ϕ(x, y, z, t), where u = ∇ϕ is the velocity field.Considering the slight compressibility of the sea, the velocity potential obeys the standard three dimensional wave equation, where c is the speed of sound in water.For the boundary condition at the free surface, we make use of results obtained from the detailed derivation given in Mei et al. (2009) section (1.1.2).We assume that (a) the wavelengths are long enough so that surface tension effects are negligible, (b) the atmospheric pressure is constant and (c) non-linear terms can be neglected over the spatial region of interest.Then the linearized, combined dynamic and kinematic boundary condition is Following Mei & Kadri (2018), the fault's ground motion is confined to a rigid, slender rectangular stripe of width 2b and length 2L, with a slenderness parameter ε = b/L ≪ 1 (see Figure 1), such that (2.4) In (2.4), w(x, y) defines the spatial extent of the rupture, τ defines the time the rupture is active and W 0 is the uplift velocity.To study the long distance propagation of acousticgravity waves we introduce re-scaled co-ordinates (see Mei & Kadri (2018)) (2.5) Letting ϕ = ϕ 0 (x, X, Y, z, t) + ε 2 ϕ 2 (x, X, Y, z, t) + • • • , the potential reduces to the twodimensional wave equation to leading order, with boundary conditions given by, The envelope of the radiated waves is governed by with the boundary conditions, Note that the boundary conditions for ϕ 0 and ϕ 2 at z = h here are different to those in the "no-gravity" case (ϕ 0 = ϕ 2 = 0).kx−ωt) dtdx, with ω representing angular velocity and k the wave number, the transformed potential is found to be,

Leading order
[ω 2 cos(µh) + µg sin(µh)] (3.1) where µ 2 = (ω 2 /c 2 ) − k 2 .The poles contributing to the contour integration derive from the dispersion relation ω 2 cos(µh) + µg sin(µh) = 0 in the denominator of (3.1).The first eigenvalue µ 0 is imaginary -all the rest are real.The first wave number k 0 corresponding to a gravity wave is always real.The following n N wave numbers [k 1 , k 2 , . . ., k N ] are also real, and correspond to the acoustic-gravity waves where N = ⌊(ωh/πc) + 1/2⌋.The gravity and acoustic-gravity modes are progressive waves.The next modes, n > N with wave numbers λ n correspond to decaying, evanescent modes (?).Thus, all modes satisfy the dispersion relation where, Inverting the transformation by contour integration we obtain the velocity potential, (3.3) Then picking out the real part for the propagating modes gives, (3.4) From which the pressure and surface elevation expressions can be obtained by differentiation using, Thus the pressure terms are now given by, (3.6) with surface elevation terms given by, (3.7) The expressions for pressure and surface elevation are in agreement with Stiassnie (2010) equations (3.13) and (3.14).

Long range modulation
Considering the region far from the fault, Mei & Kadri (2018) showed that for pure acoustic modes the envelopes A n (X, Y ) vary slowly, allowing the derivation of an analytical solution of the pressure.It is anticipated that the addition of gravity, would have a similar effect where the modal envelopes of both the acoustic-gravity modes (with the correction due to gravity), and the gravity mode (with correction due to compressbility) are all governed by the Schro dinger equation, which once solved explicitly we obtain where C(z) and S(z) are Fresnel integrals and, This result is identical in structure to that of Mei & Kadri (2018), though here it is valid also for the gravity mode n = 0.With the inclusion of these results for A n (k n , X, Y ) the leading order term for the velocity potential in equation (3.3) is now valid in the ranges x O(lε −2 ) = O(Lε −1 ), y O(lε −1 ).

Stationary phase approximation
We now apply the stationary phase approximation for different gravity phase speed conditions.We rewrite equation (3.6) as (3.10)

Acoustic-gravity modes
Consider the acoustic-gravity modes only -i.e. the second term in (3.10) -and let the phase of mode n be denoted by Γ n (ω) with, A first derivative of Γ n (ω) with respect to ω gives, at the point of stationary phase ω = Ω n , and a second derivative yields Applying the known stationary phase formula, e.g.P.539 of Richards (2009), the contribution to the pressure arising from the acoustic-gravity waves becomes, where Θ An is the phase of A n .The results for the acoustic-gravity modes are consistent with the results for the (pure) acoustic modes by Mei & Kadri (2018) with the difference that the modes here have a correction due to gravity.

Gravity-acoustic mode
To obtain the stationary phase approximation for the contribution to bottom pressure arising from the surface-gravity wave -the first term in (3.10) -consider the phase term Γ 0 (ω) for the general (compressible) case as given by, where single and doubles primes denote first and second derivatives with respect to ω.Noting that differentiation with respect to ω yields, (3.17) The stationary phase approximation requires a second derivative of k 0 , (3.18) Equation (3.18) contains terms in µ ′ 0 and µ ′′ 0 , to obtain these we differentiate the general gravity dispersion relation, where, for simplicity, quantities with tilde were normalised with lengthscale h and timescale h/g.
Following Richards (2009) and with ω = Ω 0 at the point of stationary phase, the pressure contribution arising from the surface gravity wave is given by, where Θ A0 is the phase of A 0 , and µ 0 = K 0 in this case.
For the total pressure contribution from the propagating modes, we combine (3.14) with (3.21) to give Similarly the surface elevation terms become, The forms of Γ ′′ 0 (Ω 0 ), Ω 0 , and K 0 are dependent upon any assumptions made as detailed in the three cases considered below: (1) a general solution with the compressible dispersion relation (3.19); (2) an approximate high order dispersion relation; and (3) first order shallow water approximation.The latter two assume incompressiblity.Note that for brevity cases 2 and 3 are presented in non-dimensional form.

Case 1: Compressible gravity dispersion relation
Evaluation of the surface-gravity wave contribution to surface elevation requires a method of calculation for µ 0 , Ω 0 , K 0 and Γ ′′ 0 .To obtain µ 0 we differentiate the general dispersion relation (3.19) with respect to k 0 , and make use of Γ ′ 0 (ω) = 0 at stationary phase which gives 2µ 0 x t gµ 0 tanh(µ 0 h) − gµ 0 tanh(µ 0 h) + gµ 2 0 h − gµ 2 0 h tanh 2 (µ 0 h) where Equations (3.24) and (3.25) form an implicit relationship for µ 0 , which can be solved numerically.Once µ 0 has been obtained, Ω 0 and K 0 can be derived directly from, and x + 8 t2 + x2 (3.28) Case 3: Shallow water limit In addition to the assumption of compressibility (μ 0 = k0 ) we consider the case of shallow water, i.e. tanh( k0 ) = k0 , which leads to, in agreement with Stiassnie (2010) equations (4.3a,b).To reduce to the 2D case (Stiassnie (2010)), the contribution from the envelope A 0 is removed (i.e. by setting A 0 = 1).Note that there is a factor of two magnification in the amplitude as compared to Stiassnie (2010), which is believed to be due to a typographical error in Stiassnie (2010) (see full derivation in supplementary material).

Validation
For validation purposes we use input parameters similar to those found in Mei & Kadri (2018) and Stiassnie (2010), and that are listed in Table 1.The number of acoustic modes is set at N = 10.This was shown to be a "reasonable choice" in Stiassnie (2010) and ?.The uplift velocity of 0.1 ms −1 , along with rupture duration of 10 s, implies a fault displacement of 1 m.Aside from comparison with Stiassnie (2010), further justification for using a duration of the order of tens of seconds can be found in Abdolali et al. (2017) and Grilli et al. (2013).
The current model is first validated against the theoretical solution for infinitely distributed fault, proposed by Stiassnie (2010) and a 3D numerical solver (Sammarco et al. 2013;Abdolali et al. 2015b).The later solves (2.1) with (2.2) at the surface, the movable bottom, representing the vertical uplift (2.3) and (2.4), and an outgoing Sommerfeld boundary condition at the end of numerical domain.The undistributed fault length assumption for the first set of validations allows us to use the 3D numerical solver on a vertical transect, which is computationally affordable.It allows the presence of the surface gravity waves and all available acoustic modes interacting with each other.The only constraint is the minimum grid resolution and time stepping, required for resolving the range of dominant frequencies.In this simulation, proper values are used to ensure the first 10 modes exist in the domain.
For the second set of validations, a single finite fault case is considered over a large 2D domain.The validation is conducted between the current model and a 2D numerical solver based on the Mild Slope Equation (M.S.E.) for Weakly Compressible Fluid, rigid bottom (Sammarco et al. 2013).
In the third case, a real-world multifault scenario is considered over a large domain, where the theory and the depth integrated model are compared to prove the accuracy of the theory.Then simulation on a variable bed condition is conducted to highlight the missing processes (i.e.refraction and reflection) due to the presence of seamounts and trenches.
Note that the solution, proposed by Stiassnie ( 2010) is fast, but has constant depth limit with an infinitely long fault assumption (2D model).The current theory provides a solution for the fault longitudinal extent to be finite and multi-fault condition.On the other hand, the 3D numerical solver is computationally expensive, but can take into account the entire problem without any assumption (i.e.variable depth profile).It is manageable to run it on transects, but requires massive computational resources for large 3D domains with the necessary resolution to resolve the acoustic-gravity wave field.The depth integrated model is faster than the 3D solver, but still much slower than the theoretical solution.In addition, it requires the forcing to be decomposed and solved on spectral bands.Here the validations are performed to prove the accuracy of each of the aforementioned models and theories, with proper overlaps.In other words, a coherent chain of cross-validation is performed to highlight the advantages of each method and the differences if assumptions are considered.

Bottom pressure
Consider a hydrophone station located at 1000 km along the positive x-axis.With the speed of sound fixed at 1500 ms −1 , the arrival time of the acoustic-gravity wave is approximately 670 s after the rupture.The tsunami arrives later at around 5000 s.  with proper boundary conditions at the surface and end of the numerical domain and movable bottom (bottom trace).

Surface elevation
With the inclusion of gravitational effects into the current model it is now possible to obtain surface elevation information in addition to the bottom pressure.Thus, the surface elevation results of Figure (2b) constitute new findings for the current model.This is of consequence when considering the inverse problem, since it enables evaluation of the tsunami alongside the acoustic modes, thereby reducing computation time.A remarkable correction of the tsunami amplitude is obtained (black curve) by deserting the shallow water assumption suggested by Stiassnie (2010), and instead solving the full compressible dispersion relation for µ 0 , Ω 0 and K 0 numerically.To illustrate this improvement a comparison with a full numerical solution is presented (dashed red curve).Thus, an inclusion of compressibility in the tsunami calculations provides an important correction of the amplitude and frequency (Abdolali & Kirby 2017;Abdolali et al. 2019).It is also worth noting that an accurate gravity constant should be used.
At times approaching the critical time tc = x the solution is not valid, due to constraints arising from the limitations of the method of stationary phase and approximations made in calculating stationary points (Stiassnie 2010).In this case the numerical model predicts a tsunami of peak amplitude approximately 0.6 m arriving at the critical time tc .Unfortunately all of the analytic models have a singularity at times approaching the critical time, see eq.(3.23).However, by splitting up the fault into a few parallel stripes (say over 10), each stripe has a shift in the critical time which allows calculating the contribution of most of the fault at all times.Thus, the general compressible solution can capture the main peak at tc -see Figure (2b) which serves to further validate the linear multi-fault approach.

Theoretical solution vs. Mild Slope Equation
In the previous section, the theoretical solutions for bottom pressure (3.22) have been validated against Stiassnie (2010) where the solutions exist for an infinite fault problem (y = 0).Here, the bottom pressure, calculated by the theory is validated against numerical simulations based on the Mild Slope Equation for Weakly Compressible Fluids (Sammarco et al. 2013), not only for lying on the x axis, but also for y = 0.The results are shown in  Hamling et al. (2017) discussed a fault that occurred on 14 th November 2016 in Kaikōura New Zealand.This event was reported as a "complex multi-fault rupture" -complex in the sense that at least 12 major crustal faults and extensive uplift along much of the coastline were observed.The rupture jumped between faults located up to 15 km away from each other, and individual sub-faults showed both positive and negative displacements as well as translational slipping.

Multi-fault rupture
The theory developed in section 3 is extended here to more complex situations, where two, (or more), slender faults can be combined by linear superposition.Each fault may have its own uplift duration and velocity, as well as dimension and orientation.To take account of multiple faults relative to a reference time and location, the acoustic-gravity wave component of expression (3.22) is modified as where i indexes the faults up to a maximum of M faults and t( t) is defined as, where H is the Heaviside step function, t = 0 is the time of the first fault movement, and ∆ i is the time lag for each individual fault relative to that of the first moving fault.Consider a hydrophone located on the sea-bed (z = 0) and to the right of a cluster of faults as shown in Figure (4).Then the (x, y) location of the hydrophone in each fault's co-ordinate system is given by, (5.3)

Multi-fault examples
Hamling et al. (2017) does not contain information on multi-fault geometries and timings etc. that would facilitate a validation exercise, so to link the current model with real data refer to Grilli et al. (2007) and Abdolali et al. (2017).The first paper discusses the Sumatra 2004 tsunami, and we use this to investigate agreement between the developed theory and a numerical model for acoustic-gravity waves (constant and variable depth).Since the DART network was not available at that time, we could not reliably validate the surface wave using the Sumatra event.Satellite records of surface displacement are available for Sumatra 2004 (Gower 2005), but these vary in both time and space, as the satellite moves across the Indian ocean, and thus introduce unnecessarily delicate challenges.Instead, we opted to use the Tohoku 2011 event (Abdolali et al. 2017) as reference for the surface wave validation where reliable data via the DART network is available.The DART buoys benefit from being at fixed locations while recording their time series of surface elevations.

Tohoku 2011 -surface elevation
Abdolali et al. (2017) investigated the surface gravity and acoustic-gravity wave fields produced by the megathrust Tohoku 2011 tsunamigenic event.The surface deflections generated by this event were recorded by the DART network deployed by NOAA (National Oceanic and Atmospheric Administration).The event occurred at 14:46 local time (J.S.T.) (Abdolali et al. 2017), with the tsunami waves arriving at DART buoy 21418 located at 38.735 N, 148.655E (NOAA web-site) approximately 30 minutes later (Abdolali et al. 2017).This buoy lies at a distance of about 500 km east of the epicentre, and is a good candidate for testing the surface elevation predictions made by the current model.The parameters used in the model were derived from a variety of sources.The dimensions of the fault were obtained from Encyclopaedia Britannica (https://www.britannica.com/event/Japan-earthquake-and-tsunami-of-2011).The coordinates of DART buoy 21418 referenced to the epicentre were calculated using the Haversine formula.The depth used for the calculation was an average (constant) value derived from a Google Earth transect between the epicentre and the DART buoy location.The strike angle α was taken from Okal et al. (2014), and finally the uplift and rupture duration were estimated using Figure (3) of Abdolali et al. (2017).From this figure it can be seen that the majority of the uplift had already occurred by 90 seconds after t 0the start of the rupture.The maximum uplift was 11.35 m (Abdolali et al. 2017) (8) demonstrate agreement between theory and numerics at constant depth (see the pressure field animation in the supplementary materials).Introduction of variable depth leads to the expected discrepancies between theory and the numerical model.However, even in the variable depth case, most of the important physics can be captured using the model.In    The variable depth case (red trace) shows attenuation of the signal for points along the transect past 400 km (i.e.just after the sea-mount).This shadowing effect is also apparent in Figure ( 9) where the sea-mount is to be found at approximately x = −1500 km, y = −100 km.More generally acoustic-gravity waves propagating into shallow sea depth experience frequency filtering by the water layer (Abdolali et al. 2015a;Cecioni et al. 2015).Low order modes are associated with smaller critical depths and are therefore able to propagate further onshore (Abdolali et al. 2014;Renzi 2017).These results confirm that changing sea depth cannot be ignored when making these calculations, since it affects the timing and scale of the signals measurable at any particular point.For instance, in the placement of hydrophones the water should be of a depth so as to enable recording of a large frequency range Abdolali et al. (2015c); Cecioni et al. (2015).

Displacement function
Aside from linear bottom displacement function (Eq.2.4), the sensitivity of surface gravity and acoustic waves are investigated numerically for a half-sine (Eq.5.4) and an exponential (Eq.5.5) bottom displacement function as shown in the upper panel of Figure 11 (Hammack 1973): where H is the Heaviside step function and ζ 0 is the residual displacement.For exponential displacement, ζ e (t = T ) = 2ζ 0 /3 or T = 1.11/α.The linear and exponential displacement functions result in very similar surface elevation plots, whereas their asso-Figure 7: Dynamic pressure time series for points #1, #2, #3 and #4 (from Figure 6).The black trace is the current model (constant depth), the blue trace is a depth integrated numerical model (constant depth), the red trace is depth integrated numerical but with variable depth.ciated acoustic-gravity wave plots show a difference in amplitude, with the exponential displacement function delivering a smaller amplitude.The surface elevation predicted by the current model (general compressible) for a linear displacement function is also shown.The half-sine displacement function shows a marked difference in surface elevation amplitude when compared to either of the linear, or the exponential displacement functions -approximately 50%.The half-sine displacement function is the only one of the three to have a smooth transition in velocity at t = 0.The amplitude of the acousticgravity wave produced in this case ends up larger (by t = 1600 s) than either of the linear or exponential cases, suggesting that energy is directed towards producing a larger acoustic-gravity wave at the expense of the surface wave.

Discussion
The separation of scales between acoustic and gravity waves indeed suggests analysing each wave type separately, as reported in literature.Such separation allows a comprehensive, but simplified study compromising the accuracy only slightly.However, such compromise may lead to a two-fold negative impact on the implementation of a reliable early tsunami warning system.The first, is that reducing the uncertainties is critical in the inverse problem, which can be done, with the model without adversely affecting the calculation time.The second is that an inverse approach that employs pure acoustic theory only can initially provide properties of the fault, but then calculations of the rising tsunami need to be carried out.Our model enables simultaneous calculation of all acoustic-gravity modes, including the rising gravity mode (tsunami), thus minimising the calculation time.
The current model includes a constant water depth assumption, which has implications that cannot be ignored.While the model can estimate the tsunami in the deep water, it may not be effective in describing the propagation over varying bathymetry and the shelf break for two principal reasons.The first lies in the assumption of constant depth, and thus effective techniques that take into account changes due to topography without computing the whole 3D domain need to be developed.The second reason is the neglect of sea-floor elasticity, which turns out to be important for both tsunami and acousticgravity wave arrival times (Kadri 2019).For the tsunami, neglecting elasticity results in    overestimation of the phase speed (Abdolali et al. 2018(Abdolali et al. , 2015b;;Watada 2013;Watada et al. 2014).The effect is even more dramatic for acoustic-gravity waves as they can couple to the elastic sea-floor and travel at speeds reaching 3900 ms −1 which significantly changes their arrival time (Eyov et al. 2013).
Extensions to the existing model could be developed in many ways, one of which would be to address the assumption of constant sound speed.The work of Michele & Renzi (2020) highlighted the importance of including variable sound speed profiles into practical applications such as tsunami early warning.Another area in which the model may be improved lies in the assumption of an idealised rectangular geometry for the faults.Although Table 1 of Mei & Kadri (2018) shows the rectangular fault assumption to be valid in many cases, it is obviously not valid for all.Future work could include the addition of dip and rake angles as parameters and also provision for the fault to rise with different velocities along its length, rather than the whole fault moving together with one set velocity.
.27) Equation (3.27) contains terms in µ ′ 0 and µ ′′ 0 , these are obtained from (3.20).Case 2: Third order incompressible dispersion relationNeglecting the compressibility of the water we can set μ0 = k0 in (3.19).Consideration of the first two terms in the Taylor expansion of tanh( k0 ) results in explicit forms for Ω0 Figure (2a) compares the bottom pressure signature calculated by the current model (top trace), Stiassnie (2010) (middle trace) and a 3D numerical solver, which solves equation (2.1) Bottom pressure signals predicted by current model (top), Stiassnie (2010) (middle) and numerical model (bottom).Co-ordinates are x = 1000 km, y = 0 km.Surface elevation plots generated by current model (stationary phase approximation inclusive of compressibility), Stiassnie (2010) and the numerical model.Co-ordinates are x = 1000 km, y = 0 km.
Figure (3) for a fault with dimensions of L = 100 km, b = 10 km and rise time of 2T = 10 s, with residual displacement of ζ = 1 m.As y increases, especially for y ≫ L, the signals become smaller as expected (see the pressure field animation in the supplementary materials).

Figure 3 :
Figure 3: (Upper) indicates the fault dimensions (L = 100 km and b = 10 km), the numerical domain extent and the coordinates of the virtual point observations.The time series of bottom pressure calculated from the current model (black) and extracted from the numerical model (red).Only the first mode is considered in order to keep computation time manageable.

Figure 4 :
Figure 4: Location and orientation of a slender fault cluster relative to a hydrophoneaxes and orientation of the i th fault indicated.
elevation comparison between current model (general compressible) and data recorded by DART buoy 21418 (red trace) for Tohoku 2011 event.

Figure 5 :
Figure 5: Comparison of current model against observation during Tohoku 2011.

Figure 6 :
Figure 6: Overview of area considered for the bottom pressure map.The section to the left of the black dashed line is that used in the calculations for Figure (8).The origin of x, y co-ordinates is at the earthquake epicentre (yellow star).Fault centroids are shown by blue stars and the faults delineated by rectangles.Depth below sea-level is indicated by the colour bar with the white areas at 4 km depth.The four points used to construct the time series of Figure (7) are labelled #1, #2, #3 and #4.The transect AB is shown with dashed line.
Figures (7) and (8) present results of a comparison made between the (linear) current model and a depth integrated numerical model applied to Sumatra 2004.Details of the parameters used in the model are summarised in Table (3).Computation of the bottom pressure field is for the region to the left of the dashed line shown in Figure (6).Both the time series of Figure (7) and the pressure maps of Figure(8)  demonstrate agreement between theory and numerics at constant depth (see the pressure field animation in the supplementary materials).Introduction of variable depth leads to the expected discrepancies between theory and the numerical model.However, even in the variable depth case, most of the important physics can be captured using the model.In Figure (9) we can see the superposition of pressure signals emanating from multiple slender faults with differing orientations, resulting in areas of high pressure, and areas where the signal is weaker.The pressure contours of column three in Figure (8) and the third column of Figure (9) highlight the missing processes of refraction, diffraction and interference induced by the variable sea depth and areas of localised elevation (red coloured areas Figure (6)), with refraction dominating all modes in deep water(Renzi 2017).In Figure (6) there is a transect with a sea-mount located approximately one third of the way along AB.The depth profile for this transect is shown at the top of Figure (10).Also shown in Figure (10) are pressure signals along the transect for three different times.

Figure 8 :
Figure 8: Snapshots of bottom pressure fields at t=1260 s (top row), t=1800 s (middle row) and t=2340 s (bottom row) from the current model, Eq. (5.1) (left column), numerical model for the case of constant depth of 4 km (middle column) and numerical model for the case of variable depth (right column).The domain extent is shown in Figure 6 and the boundary forcing is imposed along the dashed line -also Figure 6.The dynamic pressure variation is indicated with reference to the colour bar where white corresponds to 0 Pa.The transect AB is shown with dashed line in each subplot.

Figure 9 :
Figure 9: Maximum absolute values of the bottom pressure (P) of the acoustic wave generated by the Sumatra 2004 event during the first 1 hr since rupture.

Figure 10 :
Figure10: (Upper) The ocean profile along section AB (as shown in Figures6 and 8).Bottom pressure anomalies along transect AB at t = 1260, 1800 and 2340 s from the current model (black), numerical model with constant depth (blue) and numerical model with variable depth (red).

Figure 11 :
Figure 11: (Upper) The time series of bottom displacement for linear, half-sine (ζ s ) and exponential (ζ e ) functions.(Middle) The time series of surface elevation (η); and (Lower) Bottom Pressure signals.Coordinates are x= 1000 km, y = 0 km.

Table 1 :
Constants and parameters used in validation of current model with gravity.

Table 2 :
Constants and parameters used in the calculation of predicted surface elevation at DART buoy 21418 for Tohoku 2011 event.

Table 3 :
Parameters used for Sumatra 2004 event -ten faults in total.Includes ζ -the vertical displacement.