Fast Determination of Tsunami Source Parameters

Source parameters of tsunami waves are an essential part of any modern tsunami warning system. Recalculation of a measured time series (wave profile obtained by a seabed-based pressure sensor) in terms of initial sea surface displacement at tsunami source is among the most (or) one of the promising approaches to be applied in a warning center. The “orthogonal decomposition”, that was proposed and studied earlier by the authors, is numerically studied here. Realistic shape of sea surface displacement and digital bathymetry of the southern part of Japan are used. To study functionality of the proposed approach, wave profiles are obtained by the sensors of DONET – Dense Oceanfloor Network system for Earthquakes and Tsunamis – pressure gauge network of Japan. We stress on the quality of tsunami source parameters reconstruction as well as on time required. As observed, just a part of the first wave period is enough for robust determination of such source parameters as amplitude and total volume of water. Results of numerical tests are summarized in tables and then discussed.


Introduction
This paper is an extension of work originally presented in OCEANS'18 MTS/IEEE Kobe/Techno-Oceans 2018, May 28-31, Kobe Japan [1]. Typical tsunami warning system includes the following parts: wave generation, propagation over the water area under study, and inundation of a dry land. Aiming to achieve real time tsunami danger prediction, we address the question of fast determination of tsunami source parameters. Having such information, there are several software packages for numerical simulation of wave propagation and inundation mapping. Here we refer the reader to MOST (Method of Splitting Tsunamiofficial tool of the NOAA UnS tsunami warning centers), see [2][3][4], and TUNAMI, see [5,6].
In any case numerical simulation of wave propagation requires the initial data, say initial sea bed displacement at tsunami source. In this paper we address the question of rather fast reconstruction of the sea surface displacement at tsunami source. As an input data we use direct measurements of the wave profile (time series), obtained at any sensor system. As is well known, a number of deep-water pressure recorders are located all around Pacific Ocean, they measure and send such data via satellite channels. Among these systems we would like to mention the USA DART (Deep-ocean Assessment and Reporting of Tsunami) buoys (https://nctr.pmel.noaa.gov/Dart/), and DONET (Dense Oceanfloor Network system for Earthquakes and Tsunamis) pressure gauge network system, which is under construction in Japan (https://www.jamstec.go.jp/donet/e/).
Our goal is to recalculate these direct measurements of the water column over a single sensor (which is nothing but a curve of the wave profile) in terms of 2D function of initial sea surface displacement at tsunami source. In general mathematical statement it is not possible to solve this problem. So, we use the so-called "Calculation in Advance" inversion strategy [7]. It has been proposed to cover the subduction zone (area of possible tsunami sources) with the system of 100x50 km rectangles (typical size of the seabed displacement zone for tsunamigenic 7.5 M earthquake). In case of Alaska-Aleutian subduction zone there are 50 of such Unit Sources (UnSs), composed in two alongshore lines. Then a typical (for this region) shape of the initial sea-bed displacement, normalized with an amplitude of 0.57 m, is used at each UnSs as tsunami wave source. Propagation of such wave from all the UnSs were numerically calculated (in advance) over the entire Pacific Ocean according to the approximation of shallow water system. As a result, a database of the calculated tsunami wave profiles, generated by the source of the same initial shape at all UnSs, is available for all grid points of the Pacific Ocean. Now, having the measured (real) wave profile at certain sensor location, one can approximate this curve with the help of linear combination of calculated wave profiles from a number of UnSs, picking up at the same point of sensor location. In other words, the problem is to determine the amplification coefficients for the UnSs such that the corresponding linear combination of UnSs (say, a Composed Source -CS) provides a good approximation of the initial sea-bed displacement.
The similar approach named SIFT -Short-term Inundation Forecasting for Tsunamis (https://nctr.pmel.noaa.gov/Pdf/SIFT _3_2_2_Overview.pdf) is used by NOAA (USA) for wave height assessment in grid points of computational domain. As well as in the method presented here, on the basis of "computations in advance" of the tsunami height distribution from the set of Unit Sources the authors [7] reconstruct the CS. Unlike the method presented here, by search of set of coefficients in linear combination of UnSs, the SIFT system uses the combinatory search of values of coefficients giving a discrepancy minimum between real and calculated tsunami wave series in several tsunami detectors. Such approach requires a lot of processor time if CS consists from more, than 5 UnSs. There is no such a shortcoming of the approach here proposed, where CS restoration time practically does not depend on quantity of UnSs involved. By the way, authors of this article also took part in development of an algorithm for source inversion in the SIFT forecasting system. This approach has been many times tested against the field observations after near field and trans-oceanic events. Amazing agreement has been observed by comparison of computed wave, generated by a calculated CS, and measurements taken at various gauge stations. Nowadays, all subduction zones around Pacific, Atlantic, and Indian Oceans were covered with such systems of UnSs unit sources (https://nctr.pmel.noaa.gov/propagationdatabase.html); synthetic marigrams at every grid point from every source have been saved in a specialized database. One of the problems is to choose the aforementioned amplification coefficients in an automated mode. The orthogonal decomposition approach to handle this was proposed in [8] and then tested in [1,9]. It happens that even a part of the first wave period is enough to determine the CSthe set of amplification coefficients. In this paper we continue to study the problem of robust determination of sea surface displacement at tsunami source by a part of the wave profile.
The paper is arranged as follows. We first describe the proposed orthogonal decomposition algorithm and provide the mathematical ground behind it. In the numerical section we briefly describe digital bathymetry in use and the parameters of realistic CS, which has a feature of the reconstructed source of historical event of 1707 [10]. Then we provide the obtained numerical results in the form of a table. Finally, numerical results are discussed.

Orthogonal decomposition algorithm 2.1. Mathematical ground of algorithm
Here we provide the brief explanation of the algorithm for reconstruction of the CS. More detailed description one can find in [8].
Let f(t) be the marigramtime series, obtained at any particular sensor (DART buoy, bottom cable sensor of DONET, or other). We will use notations ( ), ( = 1, … , ) for simulated wave profiles, calculated (according to the linear or nonlinear shallow water approximation) at the same point of sensor location, by direct numerical modeling of the tsunami wave propagation. It is understood, that the subindex k refers to a tsunami wave initiated by the normalized sea surface initial profile of a given shape from the k-th UnSs. We assume that all these functions, ( ), are linearly independent. We are looking for coefficients, , for the linear combination of these functions, which provide the best approximation of the function f(t) in 2 norm: So, we consider the problem of optimal approximation of a given function, f(t), by the linear combination of the finite subset of functions from the system { ( )}. Suppose that the system is orthogonal and normalized according to: Here t0 stands for the time moment when the wave approaches the closest receiver, T is the time of the latest (of all receivers) full period on marigram which is concerned as the end of the second phase of tsunami. For sensors which are more distant off the coast than a source, the second phase of tsunami is positive, and for those which are closer to the shore, it is negative.
As is well known in the Fourier series theory, the coefficients of such optimal approximation are nothing but the Fourier coefficients of expansion of f(t) in a series with respect to { ( )} (see, for example, [11]).

Algorithm description
Based on the statement above, the following algorithm was proposed in [8] and tested in [1,9] for searching the coefficients of the linear combination allowing optimal approximation of the given function. It includes three steps. Let us transform the given system {f k } to the corresponding ortho-normal system, { } , with respect to the scalar product (2). Let Suppose that pairwise orthogonal nontrivial functions e 1 , e 2 , . . , e k-1 are already constructed. We will look for the function e k in the form: i.e. we obtain the next function, e k , as the linear combination of the new function, f k , and the already constructed ortho-normal ones, e 1 , e 2 , . . , e k-1 . The coefficients, l k,1 , l k,2 , . . , l k,k-1 , in formula (5) should be found from the conditions that the new function, e k = f k + l k,1 e 1 +. . +l k,k-1 e k-1 , is orthogonal with all the previous ones, e 1 , e 2 , . . , e k-1 , and normalized. Thus, the coefficients in (5) are calculated according to: As is well known in calculus, the optimal 2 approximation, according to criteria (1), of a given function f by the linear combination of orthogonal and normalized functions, is achieved at the Fourier coefficients, a i , of the function f(t), decomposition with respect to the system : So, to solve the original problem (1), we only need to provide the coefficients (8)  , if ≠ , ≠ + 1.
So, the proposed algorithm to approximate the given marigram, f(t), by the linear combination of synthetic marigrams, { ( )}, k=1,…,N, in order to achieve minimal 2 difference (1), is as follows.
Then the recorded marigram, f(t), is decomposed to the part of Fourier series with respect to obtained orthogonal and normalized (see (1), (2)) basis, ( ), by the formula (8).
Finally, the obtained Fourier coefficients, , are recalculated in terms of coefficients, , of the original linear combination (1), by means of (9).

Numerical experiments
Numerical experiments were arranged in the area around the south-east coast of Japan (eastern Kyushu, Shikoku islands and Kii peninsula. Geography, bathymetry and location of UnSs and CSs in this region are shown in Figure 1.

Digital bathymetry
Numerical experiments were arranged at the gridded bathymetry around Eastern Kyushu, Shikoku Island and Kii Peninsula (southern part of Japan), built on the basis of Japanese bathymetric data produced by the Japan Oceanographic Data Center (JODC), presented in Figure 1. This digital gridded bathymetry has the following characteristics: Size of array is 3000 x 2496 nodes; Grid steps are 0.003 and 0.002 degrees (which means 280.6 and 223 meters, respectively); Domain covers the area approximately between 131° E and 140° E, 30° N and 35° N.

Composition of tsunami source
In the course of numerical calculations model tsunami sources were formed of UnSs of 100х50 km in size located in two rows along the coast (Figure 1). Each of the six composed sources (CS), having the sizes of 200х100 km, is presented as the linear combination of four UnSs with given amplification coefficients (1.4; 1; 1.5; 1.1). For this set of coefficients, the maximum surface elevation in the source area is equal approximately to 3 m, and minimum water level there is estimated as -1.5 m (Figure 2). Bigger coefficients correspond to UnSs more distantly located off the coast. Numbering of CSs begins with UnSs which are at the western edge the source area ( Figure 1). Two of them (the first and the last) are indicated in Figure 1. The ocean bottom displacement in the CS used for computations is shown in Figure 2. We assume precision of 10% as the value of reasonable accuracy of the tsunami height estimation. This value allows to determine parameters of the wave near the shore with the same precision, which is enough for the rapid tsunami hazard evaluation. It means that maximum difference between the initial and calculated surface displacement is less than 10 cm for a typical initial displacement height of 1 m. Accordingly, the deviation of 1 m is considered acceptable in case the wave height is 10 m.
As a rule, a part of the wave profile from the beginning to the zero level after the highest amplitude point (see Figure 2) is sufficient to reconstruct the amplification coefficients with a particular precision. So, we tried to determine how smaller could be a part of wave profile to keep the error in reconstructed ISSD amplitude within 0.1 discretization. The location of the model tsunami source centers for numerical calculations are chosen according to estimation of the source of the major tsunami, expected, according to [10], in the near future in the Nankai subduction zone (Figure 3 is taken from [10]). The size of the model sources corresponds to a tsunamigenic earthquake of magnitude 8.0.

Positions of sources and sensor network
The virtual sensors network consists of one row of virtual deepwater wave detectors, installed along the Nankai Trough (R1-R8, indicated by green crosses in Figure 1). In addition, 3 sensors of DONET were also included to observation network system for source restoration in our numerical experiments. These sensors are indicated in Figure 4 as black squares, R9, R10 and R11. Reconstruction of a source means determination of coefficients in the linear sum of UnSs composing CS.

Numerical experiments
As a result of numerical experiments, the tsunami wave signal from each of the six locations of CS (see Figure 1) was recorded by 11 detectors, R1-R11, given in Figures 1 and 4. By means of the described method, the minimum length of a wave signal, sufficient for determination of coefficients in the linear sum of UnSs (Necessary Part of Wave Profile -NPWP) have been determined with the required accuracy. These lengths as a percentage to length of the full wave period which means complete first two phases (positive and negative) of a wave are summarized in Table 1. The first column shows receiver's number (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11). In the second column one can find the Composed Source number (1)(2)(3)(4)(5)(6). In the third and the fourth columns the length in seconds of NPWP and the full wave period are presented, respectively. Column number 5 gives the length of NPWP in percentage of the full tsunami period. The last two columns present the NPWP ending time instance and the wave arrival time to a shore. The duration of NPWP counted from the moment of wave arrival at receiver. More important for us is the end of NPWP on the record (column 6 in the Table 1). This record starts simultaneously with the initial ocean surface displacement. It is the same moment as regarded as zero in tsunami travel time. So, if NPWP (on record) is less than the tsunami travel time to the coast, then we have some time for reconstruction of CS and for estimation of the expected wave height near a shore. Numerical experiments showed that in most cases use of even a wave record only from one sensor allows to restore the source (coefficients in the linear sum of UnSs) with the required accuracy, using only a part of a wave signal (less than a half of its full period). From the results presented in Table 1 it is visible that for each CS there are several sensors, processing of a wave signal from which allows to restore the source with the required accuracy even before arrival of a wave to the nearest coast. In particular, for CS1 those are detectors R1, R2, R9, R10 and R11. For CS2 those will be sensors R1, R2, R3, R9, R10 and R11. And so on.

Discussion
The fact that among such sensors there is always a tsunami detector which is included to observation system DONET (Figures  1 and 4) is important. In these sensors (R9-R11) the first phase of tsunami on the wave record is negative (surface inundation), and then the surface elevation is coming. Finally, we can conclude that with the help of observation network described the tsunami warning system is able to estimate expected wave heights along the coast before its arrival there.
By results of the numerical experiments given in table 1 one can notice that the length of NPWP for restoration of some CSs by records of some sensors exceed 80% of the full period of a wave. Though, in most cases duration NPWP does not exceed 50 percent of the tsunami period. In particular, more than 80 percent of the full wave period is required for sufficient-quality restoration of CS having number 1 by means of sensors R1 and R6. The CS3 by means of sensors R3 and R6. And, at last, CS5 on the basis of measurements by sensor R4. For the last case the length of NWDP reaches 98.6 percent of wavelength. From Figure 1 it is visible that in all these cases the sensor is located in the direction of a long axis of CS (a rectangle 200 x 100 km) or near it. It means that such relative positioning of the CS and the tsunami detector negatively affects speed and quality of restoration of coefficients in a linear combination of UnSs composing the source. CSs are most successfully restored using wave records taken from the detectors which are located in the direction of a short axis of the source.

Conclusion
The numerical method which was used for calculating the direct initial value problem on propagation of tsunami generated by a UnSs have been many times tested on exact solutions and on real tsunami data. So, correlation between initial surface displacement at the UnSs and synthetic marigrams in receivers seems to be reliable. According to the mathematical model the increase in height of Unit Source by N times causes the proportional growth of wave height in all grid-points of computational domain.
As a result of this research it is shown that at the reasonable relative location of the area of the possible tsunami sources and deep-water tsunami detectors for correct restoration of the composed source it is enough only a part of the recorded wave signal. At the same time, a part of a wave signal, necessary for restoration, can be used by the restoration algorithm even before arrival of a wave to the coast. It gives the chance for quick solution of a direct problem of the tsunami propagation from the initial water surface displacement to estimate the expected wave height along the entire coast. Besides the carried-out numerical experiments will allow find optimal positions for installation of the new deep-ocean tsunami detectors capable to restore the composite tsunami sources more quickly.