Inversion method for initial tsunami waveform reconstruction

Introduction Conclusions References Tables Figures


Introduction
Recently, all of us were witnesses of the terrible tsunami that occurred on 11 March 2011 off the coast of Japan.The tsunami in Japan resembled the disaster in the Indian Ocean in 2004.These mega-events prompted serious efforts to address the mitigating strategies against the threats posed by tsunamis.
An increasing reliability of tsunami prediction can partially be achieved by means of numerical modeling, which makes it possible to estimate the expected propagation, as well as the run up, wave heights and arrival times of a tsunami on the coast that could be subject to risk.An important part of tsunami simulation is to gain some insight into the tsunami source.The modern tsunami early warning systems conventionally employ seismic methods to determine the source parameters.The approaches based on inversion of remote measurements of sea-level data have some advantages because seismic data are not available shortly after an event and are often imprecisely translated to tsunami data.Furthermore, tsunami wave propagation can be simulated more precisely than that of seismic waves.
Among the mathematical approaches based on inversion of near-field water-level data are the methods based on Green's functions technique (GFT) (Satake, 1987(Satake, , 1989)), least square inversion combined with the GFT (Tinti et al., 1996;Piatanesi et al., 2001) and an optimization approach (Pires and Miranda, 2001).A priori information from seismic data about a tsunami source played an essential role in the inversion method of Satake.One of the main advantages of the second group of methodologies is that they do not require a priori assumption of a fault plane solution.The first and second approaches are based on the linear shallow water theory, but the third one allows us to use the nonlinear shallow water equations or other appropriate equation sets.These methods have been widely applied in further studies of tsunami problems with various modifications (Wei et al., 2003).
Among the recent studies dealing with the inverse tsunami source problem, it is necessary to cite an approach which is successfully applied now and which is a method of direct sorting combined with the minimal residual method (Percival et al., 2011;Lavrentiev and Romanenko, 2014).In brief, this method determines the unknown coefficients of a set of the typical unit tsunamigenic sub-faults disposed at the subduction zone so that residuals between the calculating data from T. A. Voronina et al.: Inversion method for initial tsunami waveform reconstruction such combination of unit sub-faults and the real DART R data will be as small as possible in a least squares sense.
The developed numerical inversion technique based on least square inversion and truncated SVD (singular value decomposition) approach is here described to reconstruct the initial tsunami waveform (tsunami source) in a tsunami source area based on inversion of remote water-level measurements (marigrams).This inversion method was first described in its fundamentals in the Russian scientific journals (Voronina and Tcheverda, 1998;Voronina, 2004Voronina, , 2013)).Theoretical considerations of such a methodology for a linear long-wave model was discussed by Pires and Miranda (2001).
A direct problem of tsunami wave propagation is considered within the scope of the linear shallow-water theory.Numerical simulation is based on a finite difference algorithm on staggered grids.This ill-posed inverse problem of recovering initial tsunami waveforms is regularized by means of least square inversion using a truncated SVD approach.As a result of the numerical process an r-solution is obtained (Cheverda and Kostin, 2010).The proposed method allows one to control the numerical instability of the solution and to obtain an acceptable result in spite of the so-called illposedness of the problem.The efficiency of the inversion is defined by the relative errors of tsunami source reconstruction.Analysis of the singular spectrum of a matrix obtained allows one to predict the efficiency of inversion by using records produced at a given set of receivers.
One of the substantial advantages of this method is that it is completely independent of any particular source model and its distinguishing feature is the possibility to estimate the capability of a certain observation system to recover the tsunami source.
Based on the properties of the inverse operator studied numerically, we tried to answer the following questions: (1) what minimum number of marigrams should be used to reconstruct a tsunami source well enough; (2) where should the recorders (receivers) of the water-level oscillations be disposed relative to potential source area; (3) how accurately can a tsunami source be reconstructed based on recordings with a given monitoring system; (4) is it possible to improve the quality of reconstructing a tsunami source by distinguishing the most informative part of the observation system?
In order to answer these questions within the approach proposed, we have carried out a series of numerical experiments with synthetic data and different computational domains.The results of the numerical simulations have shown a promising outlook of this approach.

Models
Mathematically, the problem of recovering the initial tsunami waveform in a source area is formulated as the determination of the spatial distribution of an oscillation source us- ing remote measurements on a finite set of points (hereafter called receivers).Let us consider a coordinate system x y z and direct the axis z downwards.The plane {z = 0} corresponds to an undisturbed water surface.The curvature of the Earth is neglected.Let = {(x; y): 0 ≤ x ≤ X; 0 ≤ y ≤ Y } is a rectangular domain on the plane {z = 0}.We denote as the aquatic part of with arbitrary solid boundaries and straight-line sea boundaries.
One of the examples of spatial layout of this statement with straight-line coastal boundaries is represented in Fig. 1 and will be considered in Sect.5.2.In this case, the domain coincides with the domain .The interaction between the wave and the coast is not considered in this study.Our numerical model is based on the shallow water theory.On addition, we look for a solution only in a constrained region.The source area is assumed to be known from seismological data.Let = {(x, y) : x 1 ≤ x ≤ x M ; y 1 ≤ y ≤ y N } be a tsunami source, ⊂ ⊆ .Let η(x, y, t) be the function of the water surface elevation relative to the mean sea level.This function is considered to be a solution of the linear shallow water equation: with the initial conditions: with the completely reflecting conditions on the continental coasts: Fully absorbing boundary conditions (ABC) (Enquist and Majda, 1977) of second order accuracy are implied on the open boundaries.The acceleration of gravity is denoted as g and the wave phase velocity is defined as c(x, y) = √ gh(x, y).The tsunami wave is assumed to be triggered by a sudden vertical displacement ϕ(x, y) of the sea floor in the target domain .The variable h(x, y) is the water depth relative to the mean sea level.
The set-up inversion experiments are substantially different in the function h(x, y) which varies from h(x, y) = const to a special h(x, y) = h(x) modeling the shelf zone and, finally, to the real bathymetry of the Peru subduction zone.
The observational data are water level records which are assumed to be known at a set of points M = {(x i , y i ), i = 1, . . ., P } in the domain : One can also assume that the set of points M belongs to some line γ without self-crossing in the domain that is necessary only for theoretical purposes.

Inversion method
In short, this method is as follows.Let us denote by A the linear operator of the Cauchy problem presented by Eqs. ( 1)-( 3) and trace its solution on the line γ (s).Under an appropriate assumption on the functions h(x, y), ϕ(x, y) and the line γ (s) (this assumption does not necessarily hold in the experiments), by means of the standard technique of embedding theorems it was proved (Voronina, 2004(Voronina, , 2012) ) that the operator A : L 2 ( ) → L 2 (γ (s) × (0, T )) is compact and, therefore, does not possess a bounded inverse.Thus, Eqs.(1)-( 3) are now reduced to the following equation: where η 0 (s, t) =η(x(s), y(s), t), (x(s), y(s)) ∈ γ (s), As was shown by Kaistrenko (1972), the above inverse problem has a unique solution only if the source function allows factorization in the temporal and spatial variables.
The inverse problem in question can now be formulated as the problem of solving a linear operator equation of the first kind.Its solution will be sought for in a least squares formulation.In other words, any attempt to numerically solve Eq. (5) must be followed by a certain regularization procedure.In the present study, regularization is performed by means of truncated SVD that brings about a notion of rsolution (Cheverda and Kostin, 2010).
In brief, the notion of r-solution can be described as follows: any compact operator A can be described in a Hilbert space with a singular system {s j , g j , e j }, j = 1, . . .∞, where s j ≥ 0 (s 1 ≥ s 2 ≥ . . .≥ s j ≥ . . . ) are singular values and {g j }, {e j } are the left and the right singular vectors.A e j = s j g j and s j → 0 with j → ∞.The systems {g j }, {e j } are orthogonal.A very important property of the singular vectors is that they form bases in the data and model spaces, therefore, the solution of Eq. ( 5) can be given by the expression: η 0 • g j s j e j (x, y). (6) As one can see from Eq. ( 6), the ill-posedness of the operator equation of the first kind with compact operator is due to the fact that s j → 0 with j → ∞, i.e. one can perturb the right-hand side η 0 (s, t) in such a way that its vanishing perturbation can lead to a large perturbation of the solution.
The regularization procedure based on truncated SVD leads to a notion of r-solution given by the formula An r-solution is the projection of the exact solution of Eq. ( 6) onto a linear span of r right singular vectors corresponding to the top singular values of the compact operator A. This truncated series is stable for any fixed parameter r with respect to perturbations of the right-hand side and the operator as it is (Cheverda and Kostin, 2010).The value of r is determined by the formula r = max {k : where cond = cond(A) is set by the user restriction on the conditioning number of the matrix A. It is clear that the value of r depends on the rate of decreasing of singular spectrum of the matrix A, which is tightly bounded with the parameters of the observation system and noisiness of data.Indeed, let us assume the perturbation in the right-hand side η 0 (s, t) is known and can be written in the form then the perturbed solution will be represented as One should limit the upper index in the last sum from the time when s j is far less ε j (t) to avoid the numerical instability.It is reasonable that the larger the r is, the more informative the solution will be.Note that the fitted value of r is much less than a minimum of the matrix dimensions.The analysis of the singular spectrum of the matrix A is the key aspect in the proposed methodology because it allows one to predict the recovered source function ϕ(x, y) with a certain observation system and bathymetry.

T. A. Voronina et al.: Inversion method for initial tsunami waveform reconstruction 4 Finite dimensional approximation and r-solution
In a real situation, one can numerically resolve only a finite dimensional subsystem of Eq. ( 5) with (K × L) submatrix.Convergence of the r-solution of a finite-dimensional system of linear algebraic equations to the r-solution of an operator equation was carefully investigated by Cheverda and Kostin (2010).
To solve numerically Eqs. ( 1)-( 3) we applied a finite difference approach for its equivalent first order linear system in terms of the unknown water elevation η(x, y, t) and velocity vector V : completed by the following initial conditions: and the boundary condition on the solid boundary: and absorbing conditions on the open boundaries.This problem was approximated by an explicit-implicit four-point finite difference scheme on a uniform rectangular grid which is based on the staggered grid stencil using the centraldifference approximation of spatial derivatives.As it was mentioned above, the wave run up is not considered in this study, hence we infer the coast line when the depth h(x, y) = 50 m.
In order to obtain a system of linear algebraic equations by means of the projective method, a trigonometric basis was chosen in the model space, i.e. the unknown function of water surface elevation ϕ(x, y) is approximated in the target domain (Sect.2) by a sum of spatial harmonics {ϕ mn (x, y) = sin mπ l 1 (x − x 1 ) • sin nπ l 2 (y − y 1 )} with unknown coefficients c = {c mn }: the center of the tsunami source was believed to be at a point (x c , y c ) being the central point of the domain and l 1 = (x M − x 1 ); l 2 = (y N − y 1 ).We assumed the water level oscillations η 0 (x, y, t) are known at a set of points {(x p , y p )}, p = 1, . . ., P for a finite number of time samples {t j }, j = 1, . . ., N t , i.e.
The unknown function ϕ(x, y) was sought for according to Eq. ( 13).
Now, we assume that the dimensions of the solution and the data space are equal to: In the end, a linear algebraic system for the unknown vector c of coefficients {c mn } (ordered in any one-dimensional way) from Eq. ( 13) was obtained: where A is a matrix whose columns consist of computed waveforms in every receiver for each spatial harmonic used as the source, η 0 is a vector containing the observed tsunami waveforms.Now, the matrix A is of K by L matrix.After the SVD procedure one could obtain the singular values of the matrix A, its left and right singular vectors {g i , i = 1, . . ., L}, The value r could be fitted by analyzing the singular spectrum of the matrix A.
Then the r-solution of Eq. ( 14) is represented by the sum where {α j = (η 0 ,g j ) s j } and {g j }, {e j } are the left and the right singular vectors of the matrix A, {s j } are its singular values.
Finally, the numerical simulation includes the following steps: 1. First, we obtain synthetic marigrams in all receivers by solving the forward problem with a certain function ϕ(x, y) as a source to be reconstructed.Thus, the vector η 0 in Eq. ( 14) is obtained.
2. The computed marigrams are perturbed by a background noise, i.e. a high-frequency disturbance and appropriate filtering is applied.
4. Further, the standard SVD procedure is applied.The analysis of singular spectrum of the matrix A allows one to choose the number r varying the conditioning number of the matrix A, as was explained above, and to compute the coefficients {c mn } by solving Eq. ( 14).
6. To estimate the efficiency of the inversion experiment we use a misfit parameter which is defined as a relative L 2 -error by the formula: In the text below, the misfit parameter is denoted by err %.

Numerical experiments: description and discussion
A series of calculations has been carried out by the method proposed to clarify the dependence of the efficiency of the inversion on certain characteristics of the observation system such as the number of receivers and their location, frequency or temporal range of data and the signal-to-noise ratio.We postulate that source area is given, as it is in real cases.Synthetic data for the numerical inversion experiments presented below were computed as a solution of Eqs. ( 10)-( 12) with respective open boundary conditions.In addition, the function ϕ(x, y) was explained by the relation: where function β(x, y) makes a semi-ellipsoid with the center at the point (x 0 , y 0 ) and the radii R x and R y on the plane z = 0.The parameter α(x) is perturbation of this semiellipsoid.When α(x) = 1 the initial tsunami waveform was assumed to be a semi-ellipsoid or a semi-sphere.

Case study h(x, y) = const
Firstly, to avoid the influence of such factors as bathymetry and data noise, we consider a formal calculation domain with open boundaries, with the depth h(x, y) = h 0 , the wave phase velocity c(x, y) = √ g h 0 = c 0 and α(x) = 1.The setting of the problem allows us to consider our problem in the spectral domain and to use the r-solution method for analytical solution (Voronina and Tcheverda, 1998;Voronina et al., 2014).The basic properties of the inverse operator were numerically studied.
In short, we have shown that the inversion results are better when the aperture length is rising, increasing the upper limit of the frequency band makes the obligatory effect for the inversion quality and leads to decreasing the smearing of the recovered function and, hence, to increasing extreme values in the center of the source.The results of the inversion by a wide angle aperture coincide with the experiments on a wide linear aperture.It was shown, if receivers were uniformly distributed along the linear aperture, the efficiency of the inversion was improving when the number of receivers increases up to 15 but further increasing the receivers number was useless.The receivers-to-source distance has no effect on the inversion result in the case studied.If the number of receivers is fixed the efficiency of the inversion rises with the azimuthal coverage.We have shown that a more precise definition of the target domain relative to the size of a tsunami source leads to a decrease in the number of necessary receivers, perhaps as little as 5-7.In addition, the released energy as well as the maximum value of the recovered function increases.However, the total volume of water displaced due to the source in the domain is not varied.
The results obtained allowed us to outline the main details of the methodology proposed.The basic point is that analyzing a singular spectrum of the obtained matrix A enables us to make an assumption about the forthcoming efficiency of the inversion.
Let us see how that works in the context of the dependence of the efficiency of the inversion on the number of receivers and their azimuthal coverage.Our numerical experiments were conducted with the following parameters (in spectral domain): = {(x, y): −100 km ≤ x ≤ 100 km; 100 km ≤ y ≤ 200 km}, frequency w ∈ [0.001, 0.01] Hz; α(x) = 1 in Eq. ( 16), i.e. the type of the source is a semi-ellipsoid with R x = 50 km; R y = 25 km and the center was assumed at the point (x 0 , y 0 ) = (0.150) according to Fig. 2. The conditioning number of the matrix A was 10 8 in all calculations that is possible due noise-free synthetic data.Let us define an above assembly of the fault parameters as Model 5.1.Our purpose was to obtain acceptable results of the inversion using a minimum of marigrams, so, we made a computer simulation when the number of receivers ranged from one to five when they were placed on the circle aperture with the different aperture angle.Figure 2 shows a layout scheme of the source-receivers arrangement in Model 5.1.
We will use the following notations: (1) the conditioning number of the matrix A is designated as cond.(A); (2) the value 100 × max (ϕ(x, y)), while (x, y) ∈ is denoted by the symbol {100max}; (3) the misfit parameter is denoted by err %.
The common logarithms of singular values for these numerical experiments are plotted in Fig. 3.A sharp decrease in the singular values, when their number increases, is typical for all calculations in all cases of the study, due to the ill-posedness of the problem.
Comparing the singular spectra in Fig. 3 we expect the worst result of inversion when 1-2 receivers were used, while the best result is provided by using 5 receivers.Moreover, the behavior of singular spectra makes it possible to expect the improvement of inversion by increasing the aperture angle for every set of receivers.As will be clear below, increasing the number of records alone does not lead to a good inversion if there is insufficient 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 lower efficiency of inversion.
Indeed, in Fig. 4 one can see how the number r (the blue line) and a maximum value of the recovered function (the red line) varied on the number of receivers used (down horizontal axis) and on the aperture angle (upper horizontal axis).One can see from these graphs that increasing the number of receivers, in total, leads to a better inversion: misfit parameter (the green line) decreases by increasing the number of receivers used, as well as maximum value of the recovered function tends to a maximum value of the theoretical function.If the number of receivers is fixed, the efficiency of inversion with the azimuthal coverage that has a good matching with our previous assumption based on analyzing singular spectra.We have also shown that if the aperture angle is sufficiently wide, the influence of the conditioning number is not significant, but the inversion parameters are worse for a smaller conditioning number of the matrix A. The recovered functions for the cases discussed above are presented in Figs.5-6.

Case study h(x, y) = h(x)
How does a bottom relief and a more realistic type of the source influence the inversion result?To answer this question we have carried out numerical experiments for the model bottom topography having some basic morphological features typical of the island arc regions (see Fig. 1).
As the initial sea surface displacement, a sea floor deformation of typical tsunamigenic earthquakes with reverse dip-slip or low-angle trust mechanisms was used.Its dipolar shape was simulated according to Eq. ( 16) with the parameter α(x) = (x − x 0 + 3 • R 1 ) • (x − x 0 + R 1 /6); and R 1 = 25; R 2 = 50 (see Fig. 14, left panels).The target domain was Synthetic data for the numerical inversion experiments presented below were computed as a solution of Eqs. ( 10)-( 11).The full reflection boundary condition described in Eq. ( 12) is fulfilled on the coast line x = 0 and the absorbing boundary conditions are imposed on the open free boundaries: The function ϕ(x, y) was sought for according to Eq. ( 13) with M = 25, N = 11.Synthetic marigrams were calculated for N t = 1990 time instants.Receivers were disposed on the segment [10; 190] of the line x = 0, i.e. a maximum of the aperture length was 180.
As in the previous case with a constant depth of the calculated basin we tried to use a minimum number of marigrams in the inversion process.For this reason, the parameter P was equal to 1, 2, 3, 5.Some aspects of this study were presented in Voronina (2004).We now turn to these models to make a generalization.The point in question is the location of receivers with respect to the dipolar source.
At the initial stage, the numerical experiments were made without noisy data.Let us name the central point of our linear aperture as a midpoint.By the midpoint we mean the projection of the central point of the rectangle (or the central point of the source, which is the same) onto the aperture line (see Fig. 1).The graphs of the misfit parameter err % when one (the blue line) and two (the red line) receivers were used in terms of the position of receivers on the aperture line are plotted in Fig. 7.The red line in Fig. 7 refers to the following numerical experiments with two receivers: the first receiver was fixed at the initial point of the aperture and the second one was moving to the endpoint of the aperture seg-   ment stopping every 10 km along the line x = 0 with coordinates (0, 10 n) according to the coordinate system of Fig. 1.Hence, the length of aperture in every position is equal to 10 n, n = 1, . . ., 19.
The most important information that should be drawn from these graphs points to the presence of an obvious minimum of the functions presented.One can see points of a minimum in these graphs corresponding to the experiments with the receiver disposed at midpoint.In other words, availability of waveform from this observation point significantly improves the inversion result.The following numerical experiments confirm this inference.
A substantial decrease of the misfit parameter when receivers are placed at the midpoint can be explained by the dipolar shape of the source and, hence, the signal in this direction is mostly informative.
For this reason, in the following experiments with three and five receivers, one of receivers is always fixed at the midpoint.The conditioning number of the matrix A was equal to 100 in the experiments below.Now, synthetic waveforms were simulated as a result of the solution to the direct problem described in Eqs. ( 10)-( 11) perturbed by the background noise, consisting in a highfrequency disturbance about 5 % rate of a maximum signal amplitude over all the receivers.Admittedly, we did not obtain any appropriate result with perturbed synthetic data due the ill-posedness of the problem.However, since a tsunami wave is more lower frequency compared to the background noise, it is reasonable to apply the frequency filtration of an observation signal.In this paper, filtration is done by a method of grid function smoothing proposed by V. A. Tsetsokho and A. S. Belonosov in 1976.The description of this method can be found in Yurchenko et al. (2013).
The experiments conducted show that a higher rate of the synthetic background noise results in increasing the misfit parameter with the same conditioning number of the matrix A in spite of filtration.It is clear that increasing the conditioning number of the matrix A allows one to use a larger value r and, therefore, to obtain a more precise solution and a lesser misfit parameter.We have considered the following versions of the receivers disposition: -Case1: an observation system of three receivers.One receiver is always placed at the midpoint but two other receivers are symmetrically placed with respect to the midpoint with distances 10 n, n = 1, . . ., 9 every step.
-Case2: an observation system of five receivers.Again, one receiver is always placed at the midpoint, two pairs of receivers moved symmetrically from the center to the endpoints of the aperture, while a distance between every nearest-neighbor receivers in every pair was fixed as 40 km (the case named C2.1),20 km (the case named C2.2) and 10 km (C2.3).
There is an aperture length on the horizontal axis in Fig. 8.The parameters of the inversion for Case 1 as functions of the aperture length are presented in Fig. 8.The misfit parameter err % (the blue line) is decreasing but maximum and minimum values of the recovered function (the orange and the brown lines, respectively) converge to the corresponding values of the initial function (the red lines) at the moment when the aperture length approaches the projection of a source on the aperture line that corresponds to the yellow columns in Fig. 8.These regularities remain valid for the inversion with five waveforms.Indeed, the behavior of the 100max (the magenta line) and err % (the dark magenta line) for case C2.2 is similar to the ones of case C1 and bears witness to better results of the inversion.
Case C2.2 differs from case C2.3 by a more uniform distribution of the same number of receivers within the acting aperture segment being the projection of the source on the aperture line that leads to improving the inversion parameters.By comparing the values of parameters in the yellow and light-blue columns in Fig. 8 one can conclude that a further increase in the aperture length leads to worse results.
The importance of the receiver location at the midpoint was confirmed in these series of numerical experiments, too.As was shown, replacing the central point of a set of receivers from the midpoint results in increasing the misfit parameter and, at the same time, in decreasing the maximum value of the recovered function.In other words, we lose the most informative waveform.
As an example, the recovered functions corresponding to the inversion for Case 1 and Case 2.2 are presented in Fig. 9 (middle-right and right panels).From the graphs of singular spectra of the inversions with three and five waveforms plotted in Fig. 9 (left panel) one can expect that the inversion in the latter case will be more successful.Indeed, results of numerical experiments presented in Fig. 9 (middle, right panels) substantiate our assumption based on analyzing singular spectra.
It is now clear that the quality of inversion strongly depends on the disposition, the number of receivers and noisy data.The robust result could be obtained when, at least, one of the receivers involved is placed at the midpoint which is a projection of the major variability direction of the source.

Case study: real bathymetry of the Peru subduction zone
Finally, to illustrate some of the ideas let us consider the results obtained for the case study with the real bathymetry of the Peru subduction zone.We are interested in how distinctive features of a real bathymetry affect the inversion process.The simulation area is located from 85 to 71 • W and from 5 to 15 • S. We set up the following parameters for our calculations: the domain was the aquatic part of a rectangle = {0 ≤ x ≤ 600; 0 ≤ y ≤ 400} with piecewise-linear boundaries of the dry land, the domain { = {400 ≤ x ≤ 500; 200 ≤ y ≤ 300} (see Fig. 10).As mentioned above, the wave run up was not considered.The observed data were simulated as a solution of Eqs. ( 10)-( 12) and completed with the full absorbing boundary conditions of second order of accuracy being fulfilled on the open boundaries.
Again, as an initial sea surface elevation, the sea floor deformation of typical tsunamigenic earthquakes with reverse dip-slip or low-angle thrust mechanisms was used and it is plotted in Fig. 14 (left panel) with the following parameters: the center point (x 0 ; y 0 ) = (450; 250), maximum and minimum values of initial displacement in source area ϕ(x, y) ϕ max = 1.959 m; ϕ min = −0.67 m.As before, we tried to obtain acceptable results of the inversion using a minimum number of records.In our calculations, the function ϕ(x, y) was sought for according to Eq. ( 13) with M = 25, N = 11.These values are defined by the shape of theoretical function ϕ(x, y).Depth h(x, y) was assumed to be the real bathymetry of the Peru subduction zone.We placed 14 receivers at points in the domain which are enumerated according to Fig. 10 and are marked by the green color (•), i.e.P = 14 in Eq. ( 4).The time interval was long enough for the tsunami wave to  reach all receivers, specifically, the number of time steps was N t = 1684 in the case presented.
After specifying all the parameters we carried out the steps of the numerical simulation mentioned in Sect. 4. Synthetic marigrams were perturbed by the background noise with the disturbance rate about 3 % of a signal maximum amplitude over all the receivers.It is imperative that the filtration procedure was made after the noise pollution.Next, the standard SVD-procedure was performed and the singular spectrum of the matrix A was obtained.
As before, first of all we analyzed the singular spectrum of the matrix A. For example, the graphs of the common logarithm of the singular values of the matrix A corresponding to their numbers are presented in Fig. 11. for the inversion with three (the red line), four (the blue line), nine (the magenta line) and ten (the green line) waveforms used in the calculations.Comparing singular spectra for the inversion with some sets of three and four marigrams in Fig. 11, one can assume that such an increase in the number of receivers will lead to the deterioration of the solution.Inversion parameters corresponding to these cases are presented in Table 1 (the first and the fifth rows).
We have seen, until now, that adding new receivers always led to better solutions when a depth of the calculation basin was a synthetic function.However, this prediction is not universally true for the case with a real bathymetry.Only increasing the number of receivers without considering their azimuthal coverage could result in a worse inversion due to rising a general level of noise pollution of data.
Obviously, the parameter r should be taken only from the first interval, where the common logarithms of singular values are slightly sloping because further increases in number r leads to the solution instability.On the other hand, r should be large enough to provide a suitable spatial approximation of ϕ(x, y).From the numerical experiments it is clear that a satisfactory parameter value has r ≥ 70 (Voronina, 2012).Some results of our numerical simulations confirming these inferences are presented in Table 1.Analyzing the values of the inversion parameters in Table 1, one can see that the inversion with five marigrams (the eighth row) is better than the inversion with nine or ten marigrams (the 11th, 12th, 14th and 15th rows).To compare the inversion parameters in Tables 1 and 2 one has to look at Fig. 10.It becomes clear that the key role in improving the quality of the inversion is played by locating the monitoring system relative to the topography and source area.
Indeed, adding to the observation system the receivers with numbers {12, 13, 14} (see the 12th and the 15th row in Table 1) which were not affected by disturbance due to the reflection from the underwater vertical ledge, did not improve the solution, according to the expectations.The same is true for a set of receivers {1, 2, 3, 4} which, on the contrary, were placed in the direction where the distortion of the signal is more possible due to the features of the relief.Using receivers with numbers {5, 6, 7, 8, 9, 10, 11} one can obtain a strong improvement of the inversion.The replacement of even one or two of them leads to a solution deteri-oration (see the eighth and the ninth rows in Table 1).What are the odds?It seems reasonable to say that the matter is in the dipole shape of our source and in the existence of a certain angle between its axis and the line of the reflection (the coast line and the underwater ledge).The axis of the source in our case is directed along the axis x and is not perpendicular to the coast line like in the previous model.The latter set of receivers is distributed along the reflection ray corresponding to the direction of the strongest variability of the source.It should be remarked that the use of all the receivers does not improve the inversion because this set of marigrams involves many fewer informative waveforms that only leads to increasing the noisy pollution of the data.
The last series of the numerical experiments was aimed at finding out the most informative direction of the receivers location in terms of improving the inversion.We have chosen several sets each consisting of seven receivers differing in their location.
The upper parts of singular spectra for these numerical experiments are plotted in common logarithm scale in Fig. 12.It is also shown in Fig. 12, how one could define the number r for every settled conditioning number of the matrix A.
In Table 2 one can see the major parameters of the inversion with seven marigrams by Models V1-V4 which differ in the receivers location and the conditioning number of the matrix A.
The recovered functions for models V1.1, V2.1, V3.1, V1.2, V2.2, V3.2 are plotted in Fig. 13.It is clear, the use of the large conditioning number of the matrix A and, as a consequence, a large value of r in the liable interval leads to  , 3, 5, 7, 9, 11, 13 Maximum and minimum values of the recovered function are denoted by φmax and φmin respectively (scaled in meters).The dimension of the subspace where the exact solution was projected is symbolized by r.The conditioning number of the matrix A is denoted by cond and the misfit parameter is denoted by err %.After the inversion by model V3.2 was completed, we again solved the direct problem with the recovered and smoothed function φ(x, y) and calculated marigrams at the  same 14 points.Marigrams computed with the recovered tsunami source have a good matching with the initial synthetic ones not only in seven receivers used in the inversion process but also in all 14 receivers as well (see Fig. 15).This fact can be treated as a verification of the inversion algorithm.In addition, the coincidence of marigrams has a great importance for predicting water elevation in the area relying on data provided by some monitoring system.

Conclusions
We have applied an approach based on SVD and r-solution methods to recovering the initial tsunami waveform in a tsunami source area.
Searching for the initial water surface displacement as a series of spatial harmonics is a common practice; however, the method proposed allows one to control the numerical instability by tools of the r-solution method.The value r is essentially determined by a real bathymetry, spatial location of the observation system, the data noisiness and still it is significantly less than the dimension of the matrix obtained in the calculations.By analyzing the singular spectrum of the matrix it is possible to make a preliminary prediction of the efficiency of the inversion with a given set of the recording stations.Hence, we could make a well-aimed precomputation with varying locations of candidate receivers to obtain the best inversion possibly in a real process.This potential of the method should be kept in mind when designing a monitoring system for tsunamis.
By the numerical simulations we have shown that in order to obtain a reasonable quality of the source restoration we need at least five to seven receivers in the simulation with a real bathymetry.Reconstructing the shape of a source is improved when not only the number of receivers increases but their azimuthal coverage with respect to a source area improves as well.The receivers-to-source distance does not essentially influence the inversion result.One of the advantages of the method proposed is that it is completely independent of any particular shape of a source, however, a limitation of this method is that it is suitable only for the linear theory.
The location of receivers on direct and reflected rays corresponding to the direction of the strongest variability of the dipole source have the greatest effect for the inversion result.As is shown in the numerical experiments, it is this set of receivers that provides the best results of the inversion process.

Figure 1 .
Figure 1.The model bottom topography having some basic morphological features typical of the island arc regions.

Figure 2 .
Figure 2. A layout scheme of the source-receivers arrangement in Model 5.1.

Figure 3 .
Figure 3.The dependence of singular values of the matrix A ( in logarithmic scale) on their numbers; ρ denotes the number of receivers used in the inversion and arranged on the aperture with angle equal to α.

Figure 4 .
Figure 4.The parameters of the inversion using 1, 2, 3, 5 receivers: 100 max symbolizes maximum value of the recovered function multiplied by 100 (the red line), the values of number r (the blue line) and the misfit parameter err % (the green line).

Figure 7 .
Figure 7.The misfit parameter err % in terms of the position f receivers on the aperture: with one receiver used (the blue line); with two receivers (the red line).

Figure 8 .
Figure 8.The inversion parameters for the Case1 (three receivers) and Case2.2 (five receivers).Aperture located symmetrically with respect to midpoint.The yellow color of columns denote that the aperture is placed inside the segment of the coast line corresponding to the projection of the source on the coast.

Figure 10 .
Figure 10.Isolines of the Peru subduction zone with depth values (in km), the target domain and 14 receivers marked by the green color symbols (•).

Figure 11 .
Figure 11.Typical graphs (upper parts) of singular values in the common logarithm scale of the matrix A with respect to their numbers when the number of waveforms used in the inversion is equal to three (the red line), four (the blue line), nine (the magenta line) and ten (the green line).

Table 2 .
The inversion parameters for different sets with seven receivers.