Modeling the initial water displacement in the tsunami source area based on r-solution method is presented. This approach is independent of the earthquake parameters, because there are used only observed tsunami waveforms and a roughly estimated tsunami source area. The method proposed suppresses the negative efect of the ill-posedness of the problem determining the inevitable instability of the numerical solution. Furthermore, this approach allows one to obtain a more reasonable strategy for deploying the tsunami monitoring system. In this paper, the tsunami source of the 2013 Solomon Island tsunami event was reconstructed by the method proposed, and the deployment of DART Buoys monitoring system was examined for the eficiency to infer the tsunami source.
The mega thrust earthquakes quite often result in large tsunamis that may inflict a severe loss and pain to the population of coastal communities. Based on the experience gained from the large Indian Ocean 2004 and the Tohoku 2011 Tsunamis the humanity recognized the importance of the real time monitoring of such severe events. Since then, concerted eforts of the scientific community have been devoted to developing the numerical simulation tools and other related technologies of the tsunami waves to mitigate the adverse impacts of plausible tsunamis.
Tsunamis are very long gravity waves, with wavelengths of tens to hundreds kilometers, which exceed the ocean depth. Under such conditions, in the deep ocean their propagation can be described by the shallow-water theory. Tsunamis can be triggered by a variety of geophysical phenomena. In the first and more common case, an earthquake occurred near the sea floor, may produce a coseismic deformation that can cause a displacement in the sea floor that can, in turn, cause an initial sea surface deformation that may result in a tsunami wave. This sea surface deformation will be called an initial tsunami waveform or, simply, a tsunami source. To accurately forecast the inundation and run up in the near-field coast, where a warning should be issued no more than 20 minutes, it is necessary to gain the insight into a tsunami source at the early stages of tsunami propagation.
The tsunami waveform inversion has the advantage in determining a tsunami source, as compared with seismic waveform inversion because seismic data are often imprecisely translated into tsunami data. Furthermore, the tsunami wave propagation can be more accurately simulated than seismic waves due to the fact that bathymetry is better known than subsurface seismic velocity structure. Numerous studies deal with application of tsunami waveforms inversion for determining the tsunami source characteristics [ e.g., Satake, 1989; Tinti et al., 1996; Piatanezi et al., 2001; C.Pires and P.M.A.Miranda, 2001; Wei et al., 2003; Titov et al., 2005; Baba et al., 2009; Percival et al., 2011; Saito et al., 2010; Tsushima et al., 2012; Mulia et al., 2016].
In this paper, the inverse problem in question is treated as an ill-posed problem of the hydrodynamic inversion with tsunami waveforms, so, it imposes some restrictions on using the mathematical techniques. In other words, any attempt to solve this inverse problem numerically must be followed by a regularization procedure. To this end, the technique based on the least-squares inversion using the truncated Singular Value Decomposition (SVD) and r-solution methods ([13]) has been proposed and was first described in its fundamentals in [12], [14]. As a result of the numerical process, the so-called r-solution is a projection of the exact solution onto a linear span of the r first right singular vectors corresponding to the largest singular values of a compact operator of the direct problem. The properties of r-solution obtained are defined to a large extent by the properties of the inverse operator which were numerically investigated in [15], [16]. The direct problem of tsunami wave propagation is considered within the scope of the linear shallow-water theory. The computer simulation is based on a finite diference algorithm. The method proposed does not require any a priori information of a source, but only its general spatial localization assumed to be known from seismological data. Furthermore, this method allows one to control the instability of the numerical solution and to obtain an acceptable result in spite of the ill-posedness of the problem.
Presently, a lot of ofshore tsunami monitoring systems using submarine cabled seafloor observatory technology have been deployed in the deep ocean. The Deep-ocean Assessment and Reporting of Tsunamis (DART) buoy system, developed by the Pacific Marine Environmental Laboratory (PMEL) of the National Oceanic Atmospheric Administration (NOAA), is one of the deep-ocean tsunami observational systems [8]. The tsunami waveforms acquired by cabled ofshore ocean bottom tsunami meters are more available, free of the tide gauge response functions as well as the coastal and the harbor efects. Hence, the inversion approaches based on the deep-ocean observations can be used for a rapid estimation of a tsunami source, which, in turn, can be used as direct input for the real-time forecast of the tsunami impact.
Although a considerable attention has been given to developing the inversion methods to infer the initial tsunami waveform, a lesser number of studies has been devoted to revealing the influence of such characteristics of the monitoring system as the number and spatial distribution of the recording devices on the inversion results. In order to correlate these notions, a series of numerical experiments with synthetic data and diferent computational domains have been carried out using r-solution method ([15], [16]). As it was shown, the number r is tightly bounded with the parameters of the observational system.
The focus of this research is on the attempt of applying the regularities obtained for a more reasonable strategy of deployment of a tsunami monitoring system in reality. Results of numerical experiments are presented in the case study of Solomon Islands tsunami of 6th February 2013. 2
The tsunami wave is assumed to be triggered by a sudden vertical displacement of the sea floor. The tsunami propagation can be considered within of shallow-water theory. The tsunami source area is assumed to be known from the seismological data as a rectangle Ω , Ω ⊂ Φ ⊆ Π, where a rectangular domain Π is a calculation domain and Φ is the aquatic part of Π with the piecewise-linear solid boundaries Γ and straight-line sea boundaries. The problem is considered in an orthogonal coordinate system. The plane { z = 0} corresponds to the undisturbed water surface. The curvature of the Earth is neglected. The wave run up is not considered.
Let η(x, y, t) be a function of the water surface elevation relative to the mean sea level which is considered to be a solution of the linear shallow-water equations: completed by the following initial conditions: and the boundary condition on the solid boundary: ηt + g∇ · (hV ) = 0
V t + g∇η = 0 η| t=0 = φ (x, y),
V | t=0 = 0;
V · n = 0
as well as absorbing boundary conditions (ABC) of second order accuracy are
implied at the sea boundaries on the sides of the rectangle Π:
(
2 cηxt − ηtt + c2 ηyy| x=0 = 0; In the above equations, the vector V = (vx, vy) is the horizontal fluid velocity vector whose x- and y-components are, respectively, vx and vy, h(x, y) is the water depth relative to the mean sea level, g is the gravity acceleration, c(x, y) = √gh(x, y) is the wave phase velocity and n is the unit vector, outwardly directed, normal to the boundary, φ (x, y) is the initial water displacement defined in a tsunami source area Ω .
The inverse problem at hand is to infer the unknown initial water displacement
φ (x, y) as output while the observed tsunami waveforms as data input are
assumed to be known on a set of points R = { (xi, yi), i = 1, …, P } (below called
as receivers):
η (xi, yi, t) = η0 (xi, yi, t) , (xi, yi) ∈ R.
(
The unknown function of the water surface displacement φ (x, y) in the source area Ω was sought for as a series of spatial harmonics
M φ (x, y) = ∑
N ∑ cmn sin m=1 n=1 mπ l1 x · sin nπ l2 y for (x, y) ∈ [0, l1] × [0, l2], with unknown coeficients c = { cmn} .
In our case, the inverse problem data are the observed waveforms (marigrams) η = (η11, η12, . . . , η1Nt , η21, . . . , η2Nt , ηP 1, . . . , ηP Nt )T , ηpj = η(xp, yp, tj ) on the set of points (xp, yp), p = 1, . . . , P and at time instants tj , j = 1, . . . , Nt. Then the vector η containing the observed tsunami waveforms can be expressed as follows:
η = Ac,
where A is a matrix which columns consist of computed waveforms for every
spatial harmonic φ mn(x, y) = sin ml1π x · sin nl2π y used as initial condition to the
direct problem (
N ∑ βmjnφ mn(x, y), m=1 n=1 where ej = (β1j1, β1j2 . . . , βMjN )T . The solution obtained is stable for any fixed r with respect to perturbations of the right-hand side. The relationship between r and the singular values of matrix A as well as the conditioning number (noted as cond) of the matrix obtained by projection of the operator A in Eq. 7 onto a linear span of its r first right singular vectors can be expressed as r = max{ k : sk/s1 ≥ 1/cond} . Thus, the value of r is determined by the singular spectrum of the matrix A, and it is still significantly smaller than the dimension of the matrix obtained. A sharp decrease in the singular values, when their numbers increase, is typical of all the calculations, due to the ill-posedness of the problem. Increasing the value r leads therewith to a higher instability. On the other hand, parameter r should be large enough to provide a suitable spatial approximation of the function φ (x, y). It is clear that properties of matrix A and, consequently, the quality of the obtained solution are determined by the location and extent of the tsunamigenic area, the configuration of an observation system and the temporal extent of the signal. Some properties of the inverting operator in the context of retrieving a tsunami source were studied numerically in [15]. 4
A series of calculations have been carried out by the method proposed to clarify
the dependence of the eficiency of the inversion on certain characteristics of the
observation system such as the number of receivers and their location. The
inversion method described above was applied to the 2013 Solomon Islands event. The
6 February 2013 magnitude 8.0 Mw Santa Cruz Islands, Solomon Islands
earthquake (10.738◦S, 165.138◦E), depth 29 km, generated a tsunami that was
observed all over the Pacific region and caused deaths and damage locally. In Fig.1
the domain Π = { (x; y) : 140◦E ≤ x ≤ 185◦E; 17◦S ≤ y ≤ 13◦N } of tsunami
propagation calculated with GEBCO bathymetry (1-min resolution; available at
http://www.gebco.net/) is presented. We consider a Cartesian coordinate system
with the origin at the point (140◦E, 17◦S). Let the Ox-axis and the Oy-axis are
directed along the longitude and latitude accordingly. The tsunami source area is
a rectangle Ω = { 164.638◦E ≤ x ≤ 165.638◦E; 11.238◦S ≤ y ≤ 10.238◦S} . The
sea levels were recorded by the system of six (P = 6) DART r marked by the
white color (◦) and enumerated clockwise in Fig.1 : 1-55012; 2-55023;
3-52403;452402; 5-52406; 6-51425. The time interval was long enough for the tsunami wave
to reach all the receivers, specifically, the time step equaled 4sec, the number of
time steps was defined as Nt = 2000. In these calculations the values of
parameters M and N are empirically established as M = 15; N = 15. The matrix A is
about (225 × (2000 × p)), where p is equal to the number of tsunami waveforms
used in the inversion. Numerical simulation is based on a finite diference
algorithm and the method of staggered grids. A rectangular grid of 2700×1800 nodes
was placed over the domain Π while a rectangular grid of 61×61 nodes was placed
over the domain Ω , respectively. The epicenter of the tsunami source is assumed
to be at the node (1509, 376). The matrix A is computed with MOST(Method
of Splitting Tsunami) package [http://nctr.pmel.noaa.gov/model.html] adopted
to NVIDIA GPU ([17]). Further, standard SVD- procedure was applied to
matrix A. The analysis of singular spectrum of matrix A allows one to define the
number r and to compute the coeficients { cmn} as an r-solution of Eq.(
The series of the numerical experiments with real data were aimed to highlight the way of varying the observation system on improving inversion. One of the main factors which contribute to the dificulties is the complexity of the bottom relief with a plenty of submarine rocks and chains in the considered domain. Furthermore, it was interesting to obtain an acceptable results of the inversion using a minimum number of observed waveforms. As is known, increasing the number of receivers does not often lead to a good inversion if there is no optimal azimuthal coverage with respect to the source and, on the contrary, in real cases it turns out that the noisiness of data is raised resulting in lowering the eficiency of inversion.
First of all, the singular spectrum of matrix A was analyzed in every case. In Fig. 2, left, common logarithms of singular values of A are shown for diferent subsets of the receivers. In Fig. 2, right, zoomed plots from Fig. 2, left, are presented. Given a fixed bound on the conditioning number one can define value of r as an x-coordinate of the intersection point for the corresponding horizontal line and the singular value plot. Obviously, if the singular value plot decreases more or less smoothly up to some point, there is an opportunity to use larger r and, hence, to get the more informative solution. For the below considered receiver subsets using r > 21 appears to be impracticable due to the high level of the noisiness of the observed data which leads to the solution instability (see such example in Fig. 3 (a), (b)). −0.2 −0.4
Analysis of the singular spectra plays the key role to understand the relationship between the improvement of inversion and a change in the configuration of the observation system. It is clear that modifying in the subset of receivers r=49; φmax=2.9m; φmin=−3.6m
r=21; φmax=1.91m; φmin=−1.46m 2 10.24S results in changing the corresponding singular spectrum. Indeed in Fig. 2, the dashed line for the subset consisting of Receivers 1, 3, 5 and the dashed-dotted line for the subset consisting of Receivers 4, 5, 6 significantly difer, that is a consequence of the replacement of Receiver 1 by Receiver 4. This is in good agreement with a change in the inversion results presented in Fig. 4 (b), (c).
Analysis of the plots in Fig. 2 makes possible to expect that the worst inversion results would be obtained by using the subset consisting of Receivers 4, 5, 6 due to the more sharp decrease of its singular values in comparison to others. Indeed, comparison of the results presented in Fig. 4 confirms this assumption.
As it is shown in Fig. 2, the singular spectra of the monitoring systems
involving Receivers 1, 3, 5 (the dashed line), 1, 2, 3, 5, 6 (the dotted line) and
1, 2, 3, 4, 5, 6 (the solid line) are similar in appearance. It is possible to expect
(
The importance of azimuthal coverage with respect to the tsunami source and bathymetry features are illustrated by the inversion results for diferent receiver sets such as 1, 2, 5 (Fig. 4 (a)), 1, 3, 5 (Fig. 4 (b)) and 4, 5, 6 (Fig. 4 (c)). As is clear from the plots (a), (b), (c) in Fig. 4, the inversion results for the equipotent subsets and common the cond of the matrixes obtained are diferent. The results obtained by Receivers 1, 3, 5 are much better than for those subsets 1, 2, 5 and 4, 5, 6. The usage of Receiver 3 and Receiver 1 have significantly improved both the shape and the amplitude of the source (the plots in Fig. 4 (a) and (b) as well as in Fig. 4 (e) and (f) the plots in Fig. 4 (c) and (d)). The latter can be due to its perfect location in the direction of reflections from the submarine rock trail. On the contrary, a remote Receiver 4, surrounded by the islands, does not have any impact on the solution and, probably, only introduces additive noise. The same conclusion can be made for Receivers 6 and 2 from the comparison the results presented in Fig. 4 (a) and (e), as well as from the comparison the results presented in Fig. 4 (a) and (b). The fact is in our case, a decrease in the length of the records used in the inversion does not make any evidence of a positive efect.
The approach proposed provides a way to balance the number of the receivers and the quality of the inversion. Based on the analysis of a singular spectrum for each specific observation system one can define a maximal r which allows one to avoid the numerical instability.
Indeed, the results of numerical experiments presented in Fig.4 substantiate our assumption based on analyzing singular spectra.
Based on the carried out the numerical experiments, it is possible to conclude that the subsets including Receivers 1, 3, 5 are the most eficient to reconstruct the tsunami source by the method proposed.
After the inversion by tsunami waveforms from the Receivers 1, 2, 3, 5, 6
was completed, the direct problem was once again solved with the recovered
function φ (x, y) as initial condition (
The instability of a numerical solution of the ill-posed inverse problem in question in many instances is due to the noise in real marigrams that is a common feature in any real applications. An approach based on r-solution method allows one to control the instability of a numerical solution and to obtain an acceptable result in spite of ill-posedness of the problem. The method seems attractive from the computational point of view since the main eforts are required for calculating 15 10 ) (cm 0 e d u it lp−10 m A −25 5 ) m c ( e iltud 0 p m A −5 10 ) 5 m c ( iltdeu 0 p Am−5 −10 matrix A. If an observation system is fixed and tsunami-prone areas are defined, one can compute the matrix only once as a preliminary stage.
It is possible to make a preliminary evaluation of the eficiency of the inversion with a given set of recording stations by analyzing the singular spectrum of a relevant matrix. The results obtained allow to find the way to improve the inversion by selecting the most informative set of available recording stations. Since tsunami sources often have a dipolar shape, the location of receivers on direct and reflected rays corresponding to the direction of the strongest variability of the dipole source have the greatest efect for the inversion result. In addition, one should keep in mind that increasing the number of marigrams used for inversion does not always lead to an improved accuracy of the numerical solution. The receiver location efects the choice of number r by such a way: the better is the configuration of the observation system, the longer is a weakly decreasing part of the spectrum. Thus, the rate of the singular values descent which is most directly correlated with the receiver location should be considered as main parameter of the eficiency of the inversion.
The function recovered by the method proposed can find practical use both as an initial condition for various optimization approaches and for computer calculation of the tsunami wave propagation. It may be usefull to designing future observation systems for regions of perceived tsunami risk by providing a wellaimed precomputation with varying locations of potential sea level recorders. Acknowledgments. The author thanks Gorge Shevchenko,Ph.D., from IMGG FEB RAS for his assistance the preparation of data and Alexei Romanenko,Ph.D., from NSU for helping in calculating.