1 EMPIRICAL AND NUMERICAL MODAL ANALYSIS OF COASTAL AREAS PRONE TO TSUNAMI RISK

In this paper the Empirical Orthogonal Function (EOF) method and the k-f spectral technique are applied to extract from experimental data related to landslide-tsunamis propagating around the coast of a conical island, the wave modes that contribute to the wave field. The relevant modes are then compared with those obtained through numerical eigenanalysis of the long wave equation around the considered island. It is shown that it is possible to associate each EOF mode with some eigenvectors and the corresponding eigenfrequencies, both on the basis of the spatial shape, the wavenumber calculated along the coast and the frequency. Results confirm that landslide-generated waves propagate along the coast as trapped edge waves and the zero-th order mode is the most important.


INTRODUCTION
When tsunamis attack a coastal area, the local shape of the bottom and of the coast strongly affects the amplification of the waves and resonance phenomena can play a relevant role on the level of inundation (Bellotti et al., 2012a).These phenomena can be related to bathymetric features and to specific shapes of the coast.The more relevant effects due to the bathymetry are the continental shelf resonance and the generation of trapped edge waves.As far as coastal features are concerned, bays with regular shapes appear to be prone to the risk of harbour-like resonance phenomena.Given that the typical frequency spectrum of tsunamis is broad banded, multiple scale (i.e.regional and local) bathymetric and coastal features can have importance (Cortés et al., 2017).
In this paper, different approaches for the evaluation of the modal response of coastal areas are applied to a small island reproduced in a large-scale laboratory experiment (see Di Risio et al., 2009, Romano et al., 2013, De Girolamo et al., 2014, Romano et al., 2016).The study case is aimed at evaluating the generation and the propagation of landslide tsunamis along the coast of a small volcanic island similar to Stromboli (Tinti et al., 2005, Bellotti et al., 2009).The approaches applied and compared have been presented and discussed in details in two recent papers (Romano et al, 2013;Bellotti and Romano, 2017).These are the Empirical Orthogonal Functions method (Panizzo et al., 2006, Tolkova, 2010, Tolkova and Power, 2011, Camus et al., 2013), the k-f spectral analysis technique geophysics (e.g.Capon, 1969;Gupta et al., 1990;Lacoss et al., 1969) and the numerical modal analysis (Sobey 2006, Bellotti et al., 2012b).
The paper is organized as follows.First the experimental set-up is introduced, then the three analysis methods are described.Finally, the results are compared and conclusions are given.

EXPERIMENTAL SET-UP
Details on the experimental set-up and on the results are reported in Di Risio et al. (2009).Only a brief summary of key information is reported here.The experiments have been carried out in a large wave tank (50 m long, 30 m wide, 3 m deep) at the Research and Experimentation Laboratory for Coastal Defense of Polytechnic of Bari (LIC, Italy).A truncated conical island (base diameter equal to 8.90 m, maximum height equal to 1.20 m, see Fig. 1) made up of PVC sheets (thickness 0.01 m) sustained by a rigid steel frame was placed in the center of the tank.The water depth was kept constant at 0.80 m.The island roughly represents Stromboli volcano (Southern Tyrrhenian, Italy) if scaled down with prototype-to-model scale length equal to 1:1000.The slope (α) of the island flanks was kept constant to cot α = 3 (i.e. 1 vertical, 3 horizontal) in order to reproduce the typical slope of volcanic islands where landslides are likely to occur.The reproduced slope is almost the same as the so-called "Sciara del Fuoco" where lava flowing usually occurs during Stromboli eruptions and land-slides can generate tsunami waves (e.g.Tinti et al., 2005Tinti et al., , 2006)).To ensure the sliding of the landslide model, a 0.50 m wide flat slope, carefully and smoothly connected to the curved flanks, has been built along the path of the landslide.A solid landslide model shaped as a half of an ellipsoid has been used.It is 0.8 m long, 0.4 m wide and 0.05 m thick.In order to measure the water level around the conical island wave gauges, ultrasound water level sensors and run-up gauges were employed.Some of the instruments were kept fixed in space (i.e.run-up gauges), some others were placed on a steel frame that can rotate around the island center spanning a half of the island.In the Fig. 1 the possible position of the wave gauges placed on the movable arm are reported as dots.The run-up gauges have been used to measure the instantaneous displacements of the shoreline.While many tests have been carried out to study the properties of the tsunamis, in this paper reference is made to one specific test (see Romano et al. 2016 for details).

THE EMPIRICAL ORTHOGONAL FUNCTION METHOD
This section provides a brief description of the Empirical Orthogonal Function (EOF) method and reports the results of the analysis.The reader is referred to the available literature (Hasselman et al., 1998;Huang et al., 1998;Holmes et al., 2012;Karhunen, 1946;Loeve, 1955) for a detailed description of the method.The EOF decomposition is based on the hypothesis that the considered dynamical system is a stochastic and ergodic process with no mean trend.Consider a M x N matrix of time series Θ(x,t), with M the number of measurements in time and N the number of free surface gauges.The aim of the analysis is to get a decomposition of the original matrix Θ(x,t) in N spatial modes Γn(x), each one oscillating according to its own time series α n (t) (the index n spans both the N gauges and the N empirical spatial deformations).Provided that the mean value has been subtracted from each time series, the spatial correlation matrix C(x) of the observed free surface oscillations is calculated as: (1) The eigenvalue problem of the spatial correlation matrix can be formulated as . (2) Solution of this problem provides the N λ n and the orthogonal eigenvectors Γn(x), usually called Empirical Orthogonal Functions.EOF are the most probable deformations assumed by the dynamical system during the experiment, and define orthogonal directions along which the maximum possible variance of matrix C(x) is associated.By sorting with respect to the eigenvalues the related EOF are ordered as the most probable shapes assumed by the dynamical system in the observed time series.The time evolution α n (t) of each EOF can be obtained by projecting the original matrix of measured free surface oscillations Θ(x,t) on each eigenvector: (3) The matrix of data Θ(x,t) can now be expressed as (4) The application of the EOF method to the considered experimental dataset has been carried out considering preliminarily all the wave gauges, reported as black points in Fig. 1.However, given the special interest of this research on the wave modes trapped along the coast, only those close to the shoreline have been selected for the final analysis.These fall within 2m from the coast.Given the typical frequencies of the waves, that is of the order of 0.5 Hz, the cross-shore extent of the used wave gauges appears adequate to discriminate between zero-th and higher order edge waves.Within 2 m from the shoreline the zero-th order cross-shore structure of the surface elevation goes to zero, while for edge waves mode 1 the cross-shore structure shows the maximum at the shoreline and a crest at about 1.2 m from the shore.
Among the empirical modes identified by the method those associated with the largest eigenvalues are reported in Fig. 2 (see Bellotti and Romano, 2017, for complete analysis).Each subplot has been obtained by interpolating onto a regular grid the results at the positions of the wave gauges.Some of the modes (e.g. 1, 2, 3) are clearly associated with the high amplitude waves generated by the impact of the landslide, and take into account the waves propagating toward offshore.Other modes (4, 5, 6, 8, 9 etc.) are on the contrary clearly associated with trapped waves propagating along the coast, with a modal shape resembling that of edge waves, as discussed later on in details.Also mixed modes exist, which consider both waves at the landslide impact area as well as edge waves (e.g. 7, 10 etc.)For each mode potentially identified as edge waves, further analysis has been carried out as reported in detail in Bellotti and Romano (2017).The time function α(t) is analyzed using spectral analysis, and evaluating the peak frequency.The frequency is then associated with the mode.The modal shape along the coast is extracted and spectral analysis in the domain of the wavenumber is used to calculate a peak wavenumber that can be associated to the specific EOF mode.Finally, each empirical mode is associated with a pair of wavenumber and frequency values and can be represented on the k-f plane as in the left panel of Fig. 3.The free surface variance (i.e. the wave energy) explained by each mode is also considered in the plot, where the amplitude of the square marker is scaled according to the eigenvalue size.

THE WAVENUMBER-FREQUENCY ANALYSIS k-f
The k-f has been widely used in geophysics (e.g.Capon, 1969; I.N.Gupta et al., 1990;Lacoss et al., 1969) and can be applied to any geophysical signal measured by a spatial array of sensors (one-, two-or three-dimensional).The accuracy of the results depends on the spatial structure of the array (e.g.shape and extension of the array, minimum and maximum distance between two sensors, etc.) related to the properties of the measured signals.In this paper, the frequency domain beamforming method has been used.A brief description of this method is given in the following; the reader is referred to Yoon (2005) for further details.
The purpose of the k-f is to estimate the wavenumber frequency spectral density P(k,f) (i.e. the steered response power spectrum).Given the time series acquired by a spatial array of M sensors, placed at known positions x m = [x m , y m ] (where m = 1, …, M) in a two-dimensional space, the steered response power spectrum (hereinafter SRPS), for a specific frequency f = f 0 , is given by the following relationship (8) In equation ( 8) e is a steering vector, function of the wavenumbers vector k = (k x , k y ), W is a diagonal matrix that contains the shading weights w m (set equal to 1 as suggested by Zywicki and Rix, 2005), R is the spatio-spectral correlation matrix and (⋅) H denotes the Hermitian transpose operator.The spatio-spectral correlation matrix R(f 0 ), for a given frequency, is given by the following relationship (9) where each element G ij (f 0 ) represents the cross-power spectrum between the sensors i and j (where i, j = 1, …, M).The cross-power spectrum is defined as (10) where the generic element S(xm, f 0 ) represents the element corresponding to the frequency f = f 0 of the Fourier Transform of the signal measured at the m th sensor and (⋅)* denotes its complex conjugate.
The accuracy on estimating the wavenumber-frequency spectral density depends on the spatial structure of the array, compared to the properties of the signal.It is possible to define a wavenumber resolution Δk for a linear array as (11) where D is the extension (or aperture) of the array.Of course, a high wavenumber resolution allows isolating accurately two waves with adjacent wavenumbers.Furthermore, spatial aliasing phenomena may occur.In order to avoid it the minimum spatial lag d min between two adjacent sensors has to be smaller than a half of the smallest wavelength of the measured signal.Thus, it is possible to define the largest wavenumber that can be measured without spatial aliasing.This is the Nyquist wavenumber k N , defined for linear array as follows (12) For a fixed number of sensors there is a trade-off between spatial aliasing and wavenumber resolution, that (in this case) are substantially competing objectives.Anyway, the use of linear array with irregularly spaced sensors is one means of achieving a good balance between these two sampling issues.Although the technique is strictly valid for stationary processes, it has been applied to transient signals (e.g.Capon, 1969;I.N. Gupta et al., 1990), similar to those considered in this paper.
In Romano et al. (2013) the k-f method has been applied to several sets of wave gauges, also trying to evaluate the direction of the waves.Here results of the one-dimensional method applied to the run-up gauges is reported in Fig. 3 (right panel).The contour lines and the color map in the plot represent the SRPS.To help the interpretation several theoretical wave dispersion relations have been plotted on the figure.The two thick black lines represent the ones of the zero and first order edge wave modes, according to Ursell (1952), while the thin black line represents the 0th-order edge wave mode, obtained by Smith and Sprinks (1975), when a curvilinear shoreline is considered.

THE NUMERICAL MODAL ANALYSIS
Numerical solution of the homogeneous long waves equation around the island has been used to calculate the modes of the free long waves around the island.A cartesian reference frame is adopted, with origin at the center of the island; r is used to represent the distance from the center.The computational domain is limited by two circles, one representing the shoreline (r=r s ), one representing the separation from the domain and the open sea (r=r 0 ), where r s =2.1 m, r 0 =7.0 m.
While analytical solutions of the tsunami propagation around the island exist (e.g.Tinti andVannini, 1994, 1995;Renzi and Sammarco, 2010;Renzi and Sammarco, 2016) here a direct numerical solver (see Bellotti et al., 2012a, b) that is applicable to more complicated geometries is used.The model solves the classical homogeneous long waves equation in two horizontal directions: where η is the water free surface elevation, g the gravity acceleration, t is the time, h is the water depth.Differentiation is indicated with subscripts.Solutions that are purely oscillatory in time are found (see also Sobey, 2006) as , where ω is the real modal frequency, X(x,y) is the real modal spatial structure and i is the imaginary unit.As far as the offshore boundary conditions are concerned it was previously noted (Bellotti et al., 2012a) that approximate radiation boundary conditions can be incorporated in the eigenproblem, resulting in complex modes and modal frequencies, that can represent partially propagating waves, gradually damped in time due to the radiation towards the open sea.Here however a nodal line (e.g.η=0) is applied along the boundary that separates the computational domain from the open sea.This results in real eigenmodes and eigenvalues that are therefore directly comparable with those obtained by the EOF decomposition.At the shoreline, a minimum water depth of 0.5 cm is granted.The boundary conditions are as follows: η n =0, at r=r s and η=0, at r=r 0 , respectively applied at the undisturbed shoreline and at the offshore boundary; n is the outgoing normal to boundaries.Using the proposed form of the solution, Eq. ( 13) becomes: while the boundary conditions become X n =0 at r=r s and X=0 at r=r 0 .The problem admits real solutions for ω and a set of real orthogonal eigenvectors X(x,y), representing purely standing waves.The Finite Element Method is used to convert Eq. ( 14) and the boundary conditions into a system of linear equations, while the solution of the eigenproblem is obtained by applying a standard eigensolver searching for the 50 eigenmodes with the lower frequency.
Some of the eigenmodes are reported in Fig. 4 along with the frequency as calculated from the eigenvalues.Each eigenmode has been further processed in order to calculate the wavenumber along the coast, using the same procedure applied for the EOF modes.Therefore, each eigenmode is associated with a couple of frequency and wavenumber and represented in the k-f plane in Fig. 5. On the plot the eigenmodes group along lines representing the modes.The mode n=0 is represented by the markers located at the lowest values of the frequency.They align closely to the mode n=0 theoretical curve.The mode n=1 on the contrary do not seem to closely follow the theoretical one, that however is valid for a straight coast.Also results for modes n=2, 3, 4 are considered in the eigenmodes obtained by the numerical solver.

COMPARISON AND DISCUSSION
The three methods presented so far are compared in this section.First, the comparison is carried out using the k-f plane, then also a discussion on the shape of the modes obtained by the EOF and the numerical methods is carried out.
As far as the k-f plane comparison is concerned, Fig. 6 reports the position of the modes obtained using the EOF and numerical techniques.Also the color map representing the SRPS calculated by the k-f analysis is depicted.It appears that the peak energy has a frequency of the order of 0.5 Hz and a wavenumber of about 4 rad/m, corresponding to a wavelength of the order of 1.5 m.A very good overall agreement exists between the EOF and k-f results.While the EOF provides only some modes and therefore provides a discontinuous picture of the energy distribution in the k-f plane, the k-f analysis provides the wave energy content at much more (but still discrete) values of the frequency and of the wavenumber, covering the whole area of the plot reported in the figure.However, while the former provides also the modal shape, i.e. the surface elevation in the x,y space, the latter does not calculate any parameter related to the wave shape.The results of the numerical modal analysis appear to be in quite good agreement with the EOF and k-f results.Only the n=0 mode is of interest here, thus the set of modes lying in the plot in the lower region are to be considered.In general, the frequency of the predicted modes appears to be larger than that obtained by the other methods.It is important to stress that the modal analysis calculates the possible modes that can be excited by source terms, but does not provide in this case a measure of the probability that each specific mode will be important or not in reproducing the whole wave field.The dots in the figure only represent possible wave modes.
The comparison between the EOF modes and those obtained from the numerical modal analysis is reported for some selected results in Fig. 7 (see Bellotti and Romano, 2017, for full report of the results).The figure is organized using three columns: the left one reporting the EOF eigenvectors, the middle one the corresponding numerical mode and the right one the k-f plane with highlighted the modes, as detailed in the following.The correspondence between the EOF and the numerical modes has been evaluated on the basis of the position of each of them in the k-f plane.For each numerical mode of interest, the EOF whose wavenumber and frequency are most similar to those corresponding to the numerical mode has been selected.The right panels report, using a black filled square marker the considered EOF, while the red filled circular marker indicates the position of the corresponding numerical mode.For reference gray circles and square markers represent respectively all the numerical results and EOF.In general, the shape of the EOF and numerical modes associated on the basis of the position in the k-f plane is quite similar.
According to the two methods for the data analysis (EOF and k-f), it appears that the wave energy is distributed, in the k-f plane, around the frequency dispersion relationship typical of n=0 edge waves.The shape of the more relevant EOF clearly resembles that of the lowest order edge waves, being the surface elevation large close to the coast and decreasing fast (exponentially) toward the sea.It can therefore be concluded that edge waves are very important in the landslide generated waves propagation around the island.

COMPARISON AND DISCUSSION
In this paper the EOF method and the k-f analysis have been applied to extract, from laboratory data referring to landslide-tsunamis, the wave modes that are important in wave propagation and inundation of the coast.Analysis of the EOF results has provided a representative wavenumber and frequency for each mode, thus allowing comparison with the k-f spectral analysis.It has been found that the two techniques, although different in principle, provide results in close mutual agreement.The EOF method, however, also provides the shape of each mode, thus allowing a direct physical interpretation of the nature of the wave systems.Specifically, it has been found that the waves propagate along the coast as the superposition of modes, each showing large amplitude close to the shore and decreasing quickly in the offshore direction.To further give a physical interpretation of the wave modes, numerical modal analysis of the free long waves around the considered geometry has been carried out.These numerical results have been used to build discrete frequency dispersion relationships of edge waves around the specific coast and bathymetry used in the experiments.Some of the numerical modes have been found to closely resemble the EOF modes, that have been finally identified as zero-th order edge waves.
It is widely known that the amplification of long earthquake generated tsunami along the shore, largely depends on the natural modes of the considered coastal areas.This was, among others, pointed out by Tinti and Vannini, (1995), Tolkova andPower, (2011), Yamazaki andCheung, (2011);Bellotti et al., (2012b); Cheung et al., (2013).However, the role of the trapped wave modes for shorter landslide-tsunami was not completely addressed in the past.It is worth to cite Lynett and Liu, (2005) who suggested, on the basis of numerical computations, the existence of edge waves also for landslidetsunamis.On the basis of the findings of the present paper it finally appears that edge waves play a very relevant role, especially for what concerns the propagation at long distances from the generation area.

Figure 1 .
Figure 1.Plan (left) and section (right) view of the island model.

Figure 2 .
Figure 2. The 12 EOF modes with larger eigenvalues.Horizontal and vertical axis represent x and y respectively.Red indicates positive values of the free surface elevation, blue indicates negative values.

Figure 3 .
Figure 3. Left panel: EOF modes representation in the k-f domain; the square markers are scaled according to the size of the eigenvalues.Solid lines represent frequency dispersion relationship of the theoretical edge waves along a straight plane beach for modes n=0, 1.Right panel: energy distribution in the k-f plane for waves propagating alongshore.

Figure 4 .
Figure 4. Modal shapes and frequencies as obtained from the numerical modal analysis.Horizontal and vertical axis represent x and y respectively.Red indicates positive values of the free surface elevation, blue indicates negative values.

Figure 5 .
Figure 5. Frequency and wavenumber of the eigenmodes obtained by numerical analysis plotted in the k-f plane.

Figure 6 .
Figure 6.Frequency and wavenumber of the eigenmodes obtained by numerical analysis plotted in the k-f plane.

Figure 7 .
Figure 7.Comparison between some EOF modes (left column) and numerical modal analysis modes (central column); each mode is represented in the right column plots using a red dot (modal analysis) and a filled black square marker (EOF).For reference all the results of the modal analysis are reported in each plot using gray circles, while all the EOF results are represented as black squares.