Applying the R-solution method for designing a tsunami observational system

One of the promising methods of the early warning of a tsunami is obtaining data of the wave heights based on the numerical solution of the inverse tsunami problem by using the truncated singular value decomposition method as a variant of the least-squares method. The problem is considered within the framework of the linear theory of wave propagation. The technique proposed allows one to avoid the inevitable instability of the numerical solution. It is possible to choose the most informative directions for the placement of the observation stations, which is based on the analysis of the energy transfer by the spatial modes generated by each right singular vector. As it has turned out, the best location of the stations is closely related to the directions of the most intense distribution of the tsunami energy. One of the significant advantages of the approach presented is the possibility, without additional calculations of the tsunami wave propagation from a reconstructed source, to obtain the tsunami wave heights at the points at which there are no observations but which are associated when calculating the matrix of the direct problem operator. The implementation of the approach proposed of the actual event of the Chilean Illapel Tsunami of September 16, 2015, is presented.


Introduction
Improving the reliability of tsunami forecasts can be achieved in part through numerical simulations that allows one to estimate the expected propagation and run-up, wave heights and arrival times of tsunamis in protected areas. The key issue in assessing the possible characteristics of wave manifestation on the coast remains the initial conditions, i.e. the initial form of a water displacement (hereinafter referred to as the tsunami source) in the area of a seismic source. The tsunami waveform inversion is currently a widespread approach to reconstructing a tsunami source because seismic data are often inaccurately translated into tsunami data, and the tsunami wave propagation can be simulated more accurately that of seismic waves. The main mathematical approaches to reconstruct a tsunami source by solving the inverse tsunami problem are based on the numerical representation of the Green function [1], on the conjugate method [2] and on the minimum residual method [3]. These approaches are widely used in practice.
In this paper, the problem of reconstructing an initial tsunami waveform (below referred to as the tsunami source) from measurements of the incoming waveforms in a series of remote recorders is posed as an inverse problem of the mathematical physics. The ill-posed inverse problem of reconstructing an initial tsunami waveform is regularized by means of the least-square inversion using the truncated SVD approach. The operator of the direct problem is restricted to the subspace that is a linear hull of the r first right singular vectors. As a result of the numerical process, an r-solution is obtained. By changing the value r it is possible to control the instability of the solution. The approach to solving the problem of recovering a tsunami source was proposed by T.A. Voronina and V.A. Cheverda for the model cases [4,5,6], and then for reconstructing a tsunami source by using the real data [7,8,9]. The main parameters impacting on the inversion result, in addition to the noise level, are: the observation system and bathymetry of the calculation area, the dimension of the solution subspace and a set of the model functions used to represent a tsunami source.
In [8,9] there was investigated the influence of the number and azimuthal location of the tsunami recorders on solving various settings of the model. The choice of the optimal number and location of the observation stations and the way to reduce the amount of data by using an optimization approach was discussed in [10] for the case of Nankai Trough.
This study presents a methodology to select the location of the most informative part of a certain observation system aimed at improving the inversion result in the tsunami waveforms recovery. In the course of preliminary scenario calculations with synthetic data, a method proposed allows one to design an optimal version of the tsunami observation system. An example of applying the methodology proposed is presented for the real event of the Chilean Illapel Tsunami of September 16, 2015.

The method
The mathematical formulation of the problem and the approach used are described in detail in [4,5]. It is assumed that a tsunami source is described by the vertical displacement of the ocean floor which is caused by the main shock of an earthquake, i.e. the piston tsunami model. For the entire region, from the generation zone to the line of off shore gauges, the linear shallow water approximation can be used to describe the propagation of tsunami waves. The function ( , , ) describes the level oscillations of the free surface of the ocean that is obtained as the solution of the following initial-boundary value problem: Where n is the external normal vector to the solid wall, h (x, y) is the smooth function describing the bottom topography, g is the acceleration of gravity, = , , the phase velocity of the wave is defined as ( , ) = ℎ( , ). The condition of an absolutely reflective solid wall is satisfied on the coastline S. On the so-called open sea boundaries Γ, this algorithm implements the absolutely absorbing boundary conditions of the second order of accuracy. Under these assumptions, the initial displacement of the free surface ( , , ) in the domain Ω = (0, ) × (0, ) is described by some finite function ( , ) to mean a tsunami source, which is determined from the oscillations of the free surface level of the ocean known at certain set of the points M The approximation of problem (1) -(3) is carried out based on a finite difference approach. An explicitimplicit difference scheme was used built on a four-point pattern and a uniform rectangular spaced grid. The scheme has the second order of approximation w. r. t. spatial variables and the first order w. r. t. time. The unknown function ( , ) is sought for as a finite series of spatial harmonics with unknown coefficients { }: To find these coefficients, the following system is obtained where the vector consists of all the coefficients , and the vector ⃗ consists of the mareograms at all observation points and in all time counts. For the matrix obtained, an SVD analysis is performed, and a generalized normal r-solution is constructed, which is the approximately reconstructed initial tsunami waveform where = Naturally, with increasing the number r, the information content of the resulting r-solution increases, but the stability of the r-solution obtained is getting worse. The value of r is determined by the observation system, bathymetry, and noise level in the data.
3. The optimal design for placing tsunami observation systems From the first attempts to apply the approach described in [4], for a basin of a constant depth and the type of a model source in a semi-ellipsoidal form, the question of the optimal number and azimuthal coverage of the tsunami waveforms recorders was discussed. It was found that when the noise is absent, increasing the number of tsunami waveforms recorders leads to improving the quality of the inversion. Moreover, if the number of tsunami sensors is fixed, an increase in the azimuthal coverage has a significant effect. At the same time, the distances from recorders to a tsunami source do not affect the inversion result in all subsequent numerical experiments. In [6], the numerical experiments with a model bottom simulating the coastal shelf are presented. In this case, the initial water displacement was taken as a combination of positive and negative ones presented in figure 1(a). The inversion results have shown that the data obtained by the recorders located in the direction of the greatest variability of the source had the most significant meaning in reconstructing a model source. This direction corresponds to the most intense direction of energy propagation. At the same time, an increase in the linear aperture exceeding the projection of the tsunami source onto the coastline did not result in improving the inversion.
Further numerical experiments were carried out for the actual bathymetry of the Peru coast and the above model source [6]. The model observation system involved about 14 synthetic tsunami sensors, which had a sufficient azimuth coverage. It has turned out that only the location of the observation stations in the most informative directions ensures the reconstruction of tsunami waveforms not only at the observation points but also at those whose data were not used in the inversion. In the case of nonlinear coastline, the most informative direction corresponds to the rays of incidence and reflection from the shoreline of the direction of the most significant variability of a synthetic source. A comparison of the singular spectra of matrices for different observation systems makes it possible to evaluate the efficiency of these systems. The longer the first flat portion of the spectrum graph, the greater r can be used, i.e. the greater the dimension of the solution subspace and, therefore, the more informative the solution will be. However, increasing the value of r in addition to that leads to a higher instability. On the other hand, the parameter r must be large enough for a suitable spatial approximation of the function ( , ). In this way, the location of the observation system affects the properties of a spectrum.
For example, figure 1(b) shows plots of the decimal logarithms of the singular values of the matrix A that have been obtained with various observation systems for the semi-synthetic cases presented in [6]. Of course, if the sensors have a better location, fewer of them may have a more "attractive" spectrum, as is evident from comparing the red and the blue lines in figure1b. But on the other hand, calculations for actual events have shown that the spectra graphs can be close to each other, and by the above way, it is not easy to choose the optimal observation system. Second, one should analyze the changing amplitudes of the tsunami waveforms in some first modes to determine the most meaningful ones.
As the example of the Chilean Illapel Tsunami of 2015, figure 3a shows the graphs of modes with their numbers on the right, consisting of sequentially recorded tsunami waveforms in the DART buoys numbered according to the observation system shown in figure 2. First, as evident, the amplitudes of the recording wave in every DART buoy decrease with increasing the mode number. Hence, one should use no more than the first ten modes in the case. Third, one should compare the energy transferred by all modes in all tsunami sensors and evaluate distribution in the percentage of the energy carried by each mod. After that, one should select those sets of the observational stations associated with at least 2/3 of total energy carried by the tsunami wave.
In the case considered, figure 3(b) shows the distribution of the total tsunami energy among the DART buoys involved in the observation system in the above numerical experiment and caused by the incoming waves carried by all modes. The most informative data are the data from the DART buoys 32402, 32401 and 32412. Figure 4(a) shows that the first modes associated with the maximum singular values transfer more than 80% of the energy. The formula estimates the fraction of energy by the mode ∑ , = 1, … , 225.
The calculation of the tsunami wave energy propagation presented in figure 4(b) from paper [11] confirms our assumptions about the effectiveness of the above tsunami sensors. The most significant energy is transferred to the DART buoys 32402 and 32412, and lesser to the DART   [11]; the triangles denote DART buoys. Figure 5 shows an example of tsunami source reconstruction based on data observed by the DART buoys 32401 and 32402, which are associated with more than 2/3 of total tsunami energy in the case at hand, which is evident from figure 3(b). This data information content makes it possible to obtain an acceptable result of reconstructing the initial tsunami waveform. As was shown in [9], the data from DART buoys 32402, 32412, 32401, as the data from the directions of the most intense tsunami energy propagation, makes it possible to reconstruct tsunami waveforms at the points corresponding to the DART buoys 32411 and 43413 without their participation in the inversion process. Based on have been carried out the numerical experiments, we have to claim, if the data of the DART buoy 32412 do not participate in the inversion, they do not restore with any versions of the observation system. We assume this fact is due to the incorrect data from DART buoy 32412 highlighted in [12]. An example is shown in figure 6(c) when data of the DART buoys 324101 and 32402 are used in the inversion.

Conclusion
This paper proposes a technique for selecting the most informative part of a using observation system that increases the initial tsunami waveform reconstruction accuracy. Moreover, an efficient set of tsunami sensors provides, without additional calculations, computing the wave heights at the predetermined points where there were no observations. Numerical experiments with model observation systems at the preliminary computations for predicted tsunami hazard zones can help to design the placement of deep-sea sensors. These sensors must be located in the directions of the most intense propagation of the tsunami energy.