A DEPTH-INTEGRATED EQUATION FOR LARGE SCALE MODELING OF TSUNAMI IN WEAKLY COMPRESSIBLE FLUID

Pressure waves generated by fast seabed movement in weakly compressible sea water, namely hydro-acoustic waves travel at the sound celerity in water (about 1500 m/s). These waves are precursors of the counterpart long free-surface gravity waves and contain significant information on the tsunamigenic source. Measurement of hydro-acoustic waves can therefore anticipate the tsunami arrival and significantly enhance the promptness and accuracy of tsunami early warning systems. In this paper derivation of a novel depth-integrated numerical model for reproduction of hydroacoustic waves is presented and the application of this computationally efficient model on two devastating historical tsunamis of Mediterranean Sea in real bathymetry analyzed to reveal the effect of variable bathymetry. On the basis of the model results, some hints for deep sea observatory are given.


INTRODUCTION
A sudden movement of the seabed, triggered by underwater earthquakes, compresses the water column generating long gravity waves (tsunami) and pressure waves (hydro-acoustic wave).These low-frequency waves are precursors of tsunamis, since their propagation speed, i.e. sound velocity in water, is significantly larger than that of the long gravity waves.The correct detection of hydro-acoustic waves can therefore significantly enhance the efficiency and promptness of Tsunami Early Warning Systems (TEWS).State of the art of tsunami warning procedures currently relies on seismic and sea level measurements.Nevertheless seismic networks may lead to issue warnings, that can be later cancelled by the sea level records.The near-field prediction of tsunami propagation is challenging, given the limited time to spread the alarm: it is therefore unfeasible to wait the measurement of the tsunami itself before spreading the alert.The records of the faster hydro-acoustic waves can cope with the shortening decision time of spreading the alarm (Cecioni et al., 2014b).The idea of using measurements of hydro-acoustic waves for tsunami alert purpose, dates back to the work of Ewing et al. (1950).Experimental registration of low-frequency hydroacoustic waves generated by the sea-bed motion, has been measured during the Tokachi-Oki 2003 tsunami event (Nosov et al., 2007), by the JAMSTEC (Japan Agency for Marine-earth Sciences and TEChnology) observatory.Nosov and Kolesov (2007) and Bolshakova et al. (2011) have processed the in-situ bottom pressure records in order to estimate amplitude, duration and velocity of bottom displacement.Several analytical investigations have been carried out to study the physical properties of acoustic waves generated by bottom sudden displacement, clarifying that there exists a relationship between the tsunamigenic source and the hydro-acoustic waves (Chierici et al., 2010).However, applications to real cases require detailed numerical modeling.Three-dimensional models (Nosov and Kolesov, 2007) are straightforward to use, but require unrealistic computational times when applied to large-scale geographical areas.Hence the necessity of a two-dimensional model, that can retain all the physical features, yet at the same time be the basis of an efficient prediction tool.The numerical model proposed is based on the solution of the 2DH Mild Slope Equation in Weakly Compressible fluid (MSEWC), proposed by Sammarco et al. (2013).The model is applied to simulate the hydro-acoustic wave propagation generated by two main tsunamigenic destructive historical earthquakes occurred in the Mediterranean Sea: the 365 AD Crete event and the 1693 Sicily event.The depth-integrated model has been validated through comparison with the solution of the full three-dimensional weakly compressible wave problem in real-bathymetry, along vertical (2D) sections of the sea.The comparison allowed to set up some computational parameters in order to optimize the depthintegrated model.The paper is structured as follows: Section 2 describes the mathematical formulation of weakly compressible mild slope equation.Section 3 describes the application of derived MSEWC on two historical worst case scenarios occurred in Mediterranean Sea.In Section 4 discussions and conclusions are given.

HYDRO-ACOUSTIC WAVE MODEL
Consider the problem of wave propagation in a weakly compressible inviscid fluid over a mild varying sea bottom as shown in Fig. 1, where waves are generated by a moving bottom.The seabed vertical distance from the mean water level z = 0 as a function of horizontal coordinates and time is shown by z = −h(x, y, t).The continuity and momentum within fluid can be written as: Where ρ is fluid density, → U fluid velocity vector, P pressure, ν the kinematic viscosity and ν s second viscosity.∇ and ∇ 2 are respectively the gradient and the Laplacian and g is the gravity acceleration.Assume that the fluid density (ρ) is sum of a constant density ρ 0 and a perturbation ρ 1 << ρ 0 and the fluid is barotropic: where c s = ∂P/∂ρ is sound speed in the fluid.Neglecting second order terms, in the framework of irrotational fluid, Eq. ( 1) can be written in term of fluid velocity potential ( → U= ∇φ) as: Combining the time derivative of momentum with continuity, we obtain the weakly compressible wave equation in water.The free surface, bottom and outgoing boundary conditions are: where subscripts with the independent variables denotes partial derivatives.The waves must also be outgoing at infinity (L).
Expanding in a series of orthogonal functions f n (z), the classic eigenfunctions of the constant depth homogeneous problem, but with the local depth h = h (x, y, t): where the β n 's are the roots of the dispersion relation Indeed, the hypothesis of a mild slope allows us to seek the solution in the form: The forcing bottom displacement in term of eigenfunctions f n (z) can be rewritten as: Using the orthogonality property of the f n leads to obtain n th expansion coefficient: As a result, the governing equation and boundary conditions shown in Eq. ( 4) are rewritten in term of ψ n 's and f n 's as: Multyplying the first equation of ( 11) by f m , integration in depth yields to: Applying second Green's identity to the right hand side (RHS) terms of (12) yields: The values of f n (z) at boundaries are: Substituting boundary values at free surface and bottom in first and second RHS term of (13) and by virtue of Leibniz's rule the second term of the LHS combines with the third term of the RHS, so that Eq. ( 13) becomes: After doing some math and neglecting higher order terms O |∇h| 2 , ∇ 2 h , the final form of the equation describing all the hydro-acoustic and tsunami mechanics in the horizontal x, y plane: where C n (x, y) and D n (x, y) are given by: The subscript n indicates that Eq. ( 16) is valid for the generic n th mode, the gravity one and the hydroacoustic ones.The mathematical problem is solved for a single frequency ω and corresponding forcing spectrum for narrow frequency band.The superimposition of the results of many equations as Eq. ( 16), one for each mode, will lead to complete modeling of the free surface fluid potential generated by a fast sea-bed motion.For incompressible fluid, i.e. in the limit of c s to infinite, the MSEWC reduces to the classical MSE extended to allow for bottom motion (Cecioni and Bellotti, 2010a;Cecioni and Bellotti, 2010b).Details on Eq. ( 16) can be found in the paper of Sammarco et al. (2013).

LARGE SCALE APPLICATION TO MEDITERRANEAN HISTORICAL TSUNAMIS
A first large geographical scale application of the numerical model based on Eq. ( 16) is presented in this section.The model is applied to simulate the gravity and hydro-acoustic wave propagation generated by two main destructive historical earthquakes, occurred in the Mediterranean sea: the AD 365 Crete event, and the 1693 Sicily event.Their seismic parameters are reported in Table 1. Figure 2 shows the bathymetric data, used in the numerical simulation.These are taken from the National Geophysical Data Center (NGDC) data base, ETOPO1.The computational domain is a portion of the Medditeranean Sea.Dashed black lines represent the open sea boundaries, where a condition of waves free exit is imposed.The larger domain 1 is related to AD 365 earthquake whereas for 1693 Sicily event, the numerical domain is restricted to the smaller domain 2. In Fig. 2 two points are indicated, namely CTS and CP, which correspond to the position of two submarine observatories, equipped with hydrophones.The Catania Test   validate and optimize the depth-integrated model.From the simulated hydro-acoustic wave signals it has been possible to identify some general characterization of these pressure waves depending on the generation mechanism, the source location, the bottom topography and the depth of the pressure recording point.Furthermore the numerically reproduced scenarios provide preliminary indications for the implementation of innovative TEWS based on hydro-acoustic wave measurements.Details on each earthquake-tsunami event and their modeling, are given in the following subsections.
The 365 AD scenario The 365 AD earthquake was an undersea earthquake occurred on July, 21 st 365 AD in the Eastern Mediterranean, with an assumed hypo-center located off western Crete, along a major thrust fault parallel to the western Hellenic trench.Geologists today estimate that the quake intensity was 8.5 on the Richter Scale or higher, causing widespread destruction in central and southern Greece, Libya, Egypt, Cyprus, and Sicily.In Crete nearly all towns were destroyed.This earthquake was followed by a tsunami which devastated the Mediterranean coasts, killing thousands of people and hurling ships 3 km inland.
This work considers a reconstructed scenario of this earthquake, in order to numerically reproduce the generated hydro-acoustic waves.Following the works of Shaw et al. (2008) and Tonini et al. (2011), we consider almost the whole rupture of the western Hellenic arc.Making use of the simplified fault model of Tonini et al. (2011), the seismic parameters are described in table 2. The reconstructed residual sea-bed dislocation, ζ 0 has been calculated by means of the Okada (1985) formulas and it is shown in Fig. 3.The transient bottom motion is described by: where H is the Heaviside step function, ζ 0 (x, y) is the residual displacement, τ is the duration of the sea-bed motion and the sea-bed velocity is assumed constant, as ζ 0 /τ.

The 1693 scenario
On January 11, 1693 a powerful earthquake occurs offshore the East coast of Sicily, Italy.This earthquake was preceded by a damaging fore-shock on January 9.It had an estimated magnitude of 7.2, one of the most powerful in Italian history, destroying at least 70 towns and cities and causing the death of about 60,000 people.The earthquake was followed by tsunamis that devastated the coastal villages at the Ionian Sea coasts and in the Strait of Messina (Tinti et al., 2004).The strongest effects were concentrated around Augusta, where the initial withdrawal of the water left the harbor dry, followed by a wave of at least 2.4 m height, possibly as much as 8 m, that inundated part of the town.The maximum inundation of about 1.5 km was recorded at Mascali.The identification of the tsunami source is still an open issue; Tonini et al. (2011) carried out the Worst-case Credible Tsunami Scenario Analysis on the basis of a possible earthquake and a possible submerged landslide.For the analysis of the hydro-acoustic wave generation, the reconstructed earthquake scenario assumed by Tonini et al. (2011) has been chosen.The seismic parameters of Table 3 are taken as the wave generation condition.Figure 3 shows the plan distribution of the vertical residual sea-bed motion in the overall numerical domain, which occurs downward, differently from the 365 AD Crete scenario.
The depth-integrated model is tested by comparison with the 3D solution along a vertical section of the sea.The position of section AB crossing points CTS and CP is represented in Fig. 3, while Fig. 4 shows the water depth and the vertical residual sea-bed dislocation over the section AB.From the comparison between the solution of the complete 3D equation and of the MSEWC, a good agreement is found over all the points of both sections.The results at Catania Test Site and Capo Passero on the AB section are presented.
Figure 5 represents the pressure time series (left column) and the frequency spectra (right column) at Catania Test Site submarine station (upper plots) and at Capo Passero submarine station (lower plots).The depth-integrated results in the time domain almost perfectly reproduce the solution of the 3D mathematical problem.As it can be noted by looking at the spectra, the comparison is in a good agreement for the  MSEWC reproduced frequency range (0.15 -0.3 Hz).In the depth-integrated modeling the reproduction of frequencies out of this range, and relative to the second and higher hydro-acoustic mode, increases the computational time without providing a relevant improvement in the wave simulation.
Since most of the sea-bed motion occurs at a water depth of 2 km, the generated hydro-acoustic waves oscillate at a frequency close to the cut-off value of the first hydro-acoustic mode, i.e. f 1 = 0.18Hz calculated by Eq. ( 20). .This signal is however filtered when traveling on shallower waters depth, which allow the propagation only of higher frequency waves (Abdolali et al., 2014).
Given the complex bathymetry where waves are generated and propagate, the hydro-acoustic wave spectrum presents different frequency bands where energy is localized.In order to investigate the variable depth effect on the spectrum of bottom pressure, four points along section AB have been selected at different water depths and distances from the epicenter.As can be seen in Figure 4, point A is located on the lateral part of shaking fault at a depth of 2000 m while point B is located on the earthquake zone at 2283 m of water depth, point C and D are off epicenter at 3400 and 1950 m water depth respectively.Figure 6 shows the spectra of the numerically reproduced bottom pressure.The source area covers the water depth ranging from 1800 -2300 m, the minimum frequency for the first mode is 0.16 Hz (for the depth of 2300 m).While the low-frequency wave is traveling to reach the observatory points, the sea-mounts do not let the lower frequencies approximated by Eq. ( 20) for n = 1 to pass.From the first panel of Figure 6, the observatory is located in shallower depth.As a result, the frequencies larger than first natural mode corresponding to local depth can be observed ( f 1,A = 0.19Hz).In the second panel where the observation point is located on the maximum dislocation at 2283 m of water depth, the spectral peaks are settled on the first three acoustic modes f 1,B = 0.16, f 2,B = 0.49 and f 3,B = 0.83 (dashed lines).This indicates that at a certain point in the source area, the elastic waves are formed with normal frequencies determined by the local water depth.Spreading from the source zone, the spectrum can be enriched and filtered by higher frequencies representing shallower water normal modes.To reach point C at deeper area passing no seamount higher than generation area, whole generated frequencies have been permitted to cross and reach observatory point.The spectrum is coincided with first three acoustic modes shown by solid lines corresponding to generation area depth ( f 1,S = 0.16, f 2,S = 0.49 and f 3,S = 0.83Hz).On the other hand, the first four acoustic modes corresponding to observatory depth are shown by dashed lines ( f 1,C = 0.11, f 2,C = 0.33, f 3,C = 0.55 and f 3,C = 0.77Hz).The fourth panel emphasizes that the shallower depth of obervatory (h D = 1950m has barricaded propagation of pressure waves with frequencies less than 0.19Hz ( f 1,D ).The second and third natural mode related to source area can reach the observatory point.Subsequently, the spectrum is peaked at f 2,S = 0.49 and f 3,S = 0.83Hz.In order to correctly reproduce the hydro-acoustic wave field in the numerical domains of Figure 2, a maximum element mesh size of 1.2 km has been chosen, for a total number of 2,250,000 and 1,100,000 triangular elements for AD365 and 1993 Sicily events respectively.The computations have been done for frequency range 0-0.05 Hz for surface gravity wave and 0.05 -0.5 Hz for hydro-acoustic ones.Solutions were obtained using a high-speed computer equipped with 12 cores i7 3.20 GHz CPU and 64 GB RAM.The MUMPS solver (MUltifrontal Massively Parallel sparse direct Solver) has been used.The results of this large scale simulation are shown in Figures 7 for AD 365 and Figures 8 for 1693 Sicily events in term of free-surface elevation, η(x, y, t).At the time t = 0 the earthquake occurs, then the hydro-acoustic wave travels at the sound celerity in water, 1500 m/s, while the surface gravity wave front travels with V g = gh in the Mediterranean sea.Snapshots of the free-surface elevations for gravity waves (n = 0 and hydroacoustic modes (n ≥ 1) are represented with 10 and 2 minute time interval respectively.For the case of AD 365 event, a complicated pressure wave perturbation is formed at the free surface and propagates with the sound celerity in water (panel a).It covers the entire domain 12 minutes after the event whereas the surface gravity waves reach Italian coasts 1 hr after (panel b).The hydro-acoustic wave field shows that these waves tend to propagate towards deeper areas; maximum values of hydro-acoustic wave amplitude are in deeper water close to the generation area.For 1693 Sicily scenario shown in Figure 8  for Catania and Capo passero observatories, as shown in Figure 2.Each spectrogram is normalized by dividing the maximum value among two points in order to show the ratio of wave amplitudes at different distances and depths.The local first acoustic modes f 1 are shown by horizontal lines.The arrival time of low-frequency hydro-acoustic signals, traveling with sound celerity in water from the assumed fault is correctly estimated by the model.For the case of AD 365, at both points the hydro-acoustic signal arrives around 5 minutes after the earthquake occurrence, while the long free-surface tsunami waves arrive after about 40 minutes, assuming a constant velocity of gh in a mean water depth of 2.5 km.Since most of the sea-bed motion occurs at a water depth of 3 km, the generated hydro-acoustic waves oscillate at a frequency close to the cut-off value of the first hydro-acoustic mode, i.e. f 1 = c s /4h = 0.12 Hz.This signal is however filtered when traveling on shallower water depths, which allow the propagation only of higher frequency components (Kadri and Stiassnie, 2012).The model results, shown in term of frequency spectra in Fig. 9, confirm this water depth filtering effect: the CTS is located at a water depth of 2 km and therefore the cut-off frequency for propagating wave is f 1 = 0.19 Hz; the CP site is located at 3.4 km of water depth and records wave oscillation of 0.11 Hz.On the contrary, for the case of 1693 Sicily event , the observatories are at water depth deeper than generation area and close to shaking falult where there is no sea-mount between them.The CTS point lies above the earthquake, while CP point is about 100 km distant from the epicenter: hence the hydro-acoustic signal arrives after about 1 minute from the earthquake occurrence (Fig. 10).Despite the time of arrival and a larger amplitude of the waves at CTS just after the earthquake occurrence, the two pressure signals are very similar.The generated signal has a peak frequency around 0.19 Hz, since the waves are generated at a water depth of 2 km, and at the same frequency are associated propagating modes at the deeper 3 km water depth (CP).

CONCLUSION
During history, the Mediterranean sea has experienced catastrophic submarine earthquakes and is potentially hosting active faults.Tsunamis triggered by seismic activities can therefore be destructive.Implementation of a fast Tsunami Early Warning System is mandatory and for the new generations of TEWS, precursor component of tsunami waves like hydro-acoustic waves should be taken into account.The importance of immediate and accurate alarm is magnified since the coastal areas are densely populated and the travel time of tsunamis are very short due to closeness of source and destination, i.e. the travel time of tsunami wave towards the coast is in the order of few minutes up to one hour.The proposed depthintegrated model reproduces the mechanism of generation and propagation of surface gravity and hydro-  The numerical model results highlight that the sea-bed motion energy is transferred to hydro-acoustic waves mainly concentrated at the natural acoustic modes.Close to the generation area the frequency spectra clearly show energy peaks for each acoustic mode; spreading off the generation area, the seamounts and barriers shallower than the source depth filter out frequencies larger than the cut-off values.The simulation of the 1693 earthquake scenario shows that 100 km far from the epicenter, at Capo Passero observation site, the hydro-acoustic waves recorded still allow to identify the energy associated to each acoustic mode since there is no barrier between source and observatory.The numerical simulation confirms that the first mode is the one which carries most of the energy.The hydro-acoustic modes propagate undisturbed in water layer equal or deeper than the one where they have been generated.As the waves propagate in shallower water depth, characterized by higher cut-off frequencies, the components with frequencies lower than the cut-off become evanescent: i.e. when hydro-acoustic waves propagate towards shallower sea depth, the water layer acts as a frequency filter.
To implement innovative TEWS based on measurement and analysis of hydro-acoustic signals, the complete modeling of hydro-acoustic waves in real bathymetry has proven to be extremely useful.The numerical simulations show that the hydrophones must be placed in waters deep enough to record larger frequency range and, if possible, not shielded by sea-mounts.However these instruments should be connected to the shore by submarine cables to guarantee fast data transmission, therefore their location can not be too far from the coastline: hence the model can help choosing the right positioning of the hydrophones system.For the portion of the Mediterranean sea here analyzed, the numerical simulation results suggest that offshore the Sicilian East coast, where the instruments have been placed, and offshore Greek West coast, the waters are deep enough to record and identify the hydro-acoustic waves generated by seismic activities in both the Hellenic arc and the Ionian Sea (South Italy).

Figure 2 :
Figure 2: Mediterranean Sea bathymetry.Domain 1 and 2 represent the computational domain for scenario AD 365 and 1693 respectively.Points CTS and CP are deep sea observatory stations at Catania (2000 m) and Capo Passero (3400 m)

Table 1 :Figure 4 :
Figure 4: Water depth (h) and Residual sea-bed displacement (ζ 0 ) for the vertical section AB, which crosses the two CP and CTS stations.

Figure 5 :
Figure 5: Pressure time series (left column) and their frequency spectra (right column) at CTS (a; b) and at CP (c; d).Comparison between the solution of the 3-D reference model (upper plots) and the solution of the depth-integrated equation (lower plots).

Figure 6 :
Figure 6: Sea bed pressure spectra at point A (h=1990 m), B (h=2283 m), C (h=3401 m) and D (h=1950 m) of Figure 4.The dashed lines indicate the natural acoustic frequencies corresponding to observation depth while the solid lines represent the normal acoustic modes for the source depth.

Figure 7 :
Figure 7: Snapshots of the free surface (η) hydroacoustic perturbation (a) and gravity wave (b) given by the AD 365 earthquake.t = 0 refers to the time of occurance of the earthquake

Figure 8 :
Figure 8: Snapshots of the free surface (η) hydroacoustic perturbation (a) and gravity wave (b) given by the 1693 Catania earthquake.t = 0 refers to the time of occurance of the earthquake , the hydro-acoustic waves cover the entire computational domain 8 minutes after event (panel a) and the gravity waves reaches the Greek coasts almost 40 minutes later (panel b).Spectrograms of bottom pressure are presented in Figures 9 for AD 365 and 10 for 1693 scenarios

Figure 9 :
Figure 9: The bottom pressure spectrogram at Catania Test Site point (upper panel) and at Capo Passero point (lower panel), resulting from the numerical solution of the depth-integrated Eq. (16) over the whole domain.Both plots are normalized divided by the maximum value (Scenario 365).

Figure 10 :
Figure 10: The bottom pressure spectrogram at Catania Test Site point (upper panel) and at Capo Passero point (lower panel), resulting from the numerical solution of the depth-integrated Eq. (16) over the whole domain.Both plots are normalized divided by the maximum value (Scenario 1693).
acoustic waves independently due to their different frequency ranges.Solution of the Mild-Slope Equation in Weakly Compressible fluid(Sammarco et al.;2013) is computationally efficient and accurate in comparison with the solution of fully three dimensional models.Simulation of low frequency hydro-acoustic waves in weakly compressible fluid, at real-scale domains is extremely important and lead to better understading of the relevant phenomena.The depth-integrated equation has been already validated bySammarco et al. (2013); however in this paper a first application and validation in real bathymetry has been developed.The full modeling of the hydro-acoustic waves field generated by two historical tsunamigenic earthquake scenarios in the Mediterranean Sea, has been carried out.

Table 3 :
Tonini et al. (2011)f the 1693 Sicily earthquake, as reconstructed byTonini et al. (2011), for the four segments from North to South.