A New Approach to Short-Term Tsunami Forecasting

1. Introduction

A short-term tsunami forecast and effective tsunami warning is a key problem of the tsunami service. At present, the conventional method for short-term tsunami forecasting is based on seismological information only (earthquake magnitude, time of the main shock and epicentre location). An earthquake magnitude which exceeds the established threshold value, which is different for different tsunamigenic zones, typically results in the issuing of a tsunami warning (Gusiakov, 2011). This approach, based on the “magnitude-geographical principle”, is straightforward and rather effective: at least it ensures few tsunami omissions. However, this method cannot provide sufficiently accurate estimates of tsunami wave heights, especially for specific coastal areas. Firstly, this is because the dependence of the tsunami intensity on the earthquake magnitude is far from being deterministic (Gusiakov, 2011): strong earthquakes can produce weak tsunamis (or no tsunamis at all) and, vice versa, sometimes destructive tsunamis are induced by relatively weak earthquakes. Secondly, because actual tsunami run-ups are strongly variable along the coast (c.f. MacInnes et al., 2009), the tsunami alarm can be justified for one coastal site, but false for another. Moreover, the magnitude-geographical warning method does not allow evaluation of the duration of the tsunami alarm or prediction of the arrival time of the largest tsunami wave (not necessarily the first, which is quite typical for tsunami events (c.f. Rabinovich et al., 2008)). Consequently, up to 75% of tsunami warnings turn out to be false (Gusiakov, 2011; Poplavsky et al., 1997). False tsunami warnings cause anxiety among emergency managers, government officials and the business community. These false alarms result in actual financial losses, sometimes significant, including losses due to production downtime and expenses for emergency evacuation procedures and navigation of vessels to the open sea (NOAA Tsunami, 2004).

According to the modern concept, a tsunami alarm should be declared in good time and only selectively for exposed coastal areas where the incoming tsunami will be a real threat. Each tsunami alarm should be accompanied by trustworthy information on the expected tsunami arrival time, wave heights along the coast, their frequencies and the duration for dangerous oscillations to remain in effect. However, based merely on the seismological data such information cannot be obtained. Meanwhile, emergency managers and local officials are in urgent need the effective operational tools that can provide reliable tsunami forecasts as guidance for rapid, critical decisions in which lives and property are at stake (Titov et al., 2005).

From this point of view, direct tsunami monitoring in the open-ocean by DART

DART = Deep-ocean Assessment and Reporting of Tsunami is a deep-sea tsunami monitoring system by the US National Oceanic and Atmospheric Administration (NOAA).

stations can significantly improve the situation (Titov et al., 2005; Titov,2009). These deep-ocean real-time tsunameters have been developed by the Pacific Marine Environmental Laboratory (PMEL/NOAA) for the early detection, measurement and real-time reporting of far-field tsunamis. These stations have been deployed at strategic locations throughout the Pacific Ocean near the main seismically active regions with a history of tsunami generation (Mofjeld, 2009). DART tsunami monitoring systems were expected to become the backbone of the modern Tsunami Warning System (TWS) in the Pacific Ocean.

The principal methodology of the tsunami forecast based on open ocean sea level measurements was developed during 1987-2005 (Chubarov & Shokin, 1995; Chung, et al., 1995; Korolev & Poplavsky, 1995; Poplavsky et al., 1997; Satake, 1987; Titov et al., 2005; Voronina & Tcheverda, 1998; Wei et al., 2003). However, only recent advances in tsunami measurement and numerical modelling technology have made it possible to create effective tsunami forecasting systems. The tsunami forecasting technology elaborated at PMEL is based on the integration of real-time DART measurements and state-of-the-art numerical modelling (Titov et al., 2005). This technology was applied in an operative regime for several major events during 2005-2011 and was demonstrated to be highly efficient (Japan (East Coast of Honshu) Tsunami, 2011; Tang et al., 2008; Titov, 2009; Wei et al., 2008; Yamazaki et al., 2006). The maximum wave heights and other parameters of arriving tsunamis for several US sites were predicted and were later found to be in good agreement with the actual in situ observed waves (Titov, 2009). The method described in (Tsushima et al., 2011) is similar to the above method.

An inversion method proposed in (Satake, 1987) also gave good results (Fujii & Satake, 2008, 2011), however, determination of the area of a seismic centre on aftershock area for one day after the main shock is required for computation. The latter does not allow applying them in an operating mode.

It should be noted, however, that the PMEL technology is based on a pre-computed Propagation Database to compute a quick preliminary forecast of the ocean-wide propagation of the tsunami as a linear combination of unit sources selected to represent the initial earthquake parameters. Once the actual tsunami wave reaches a DART station, inversion is performed to adjust the slip distribution of the selected unit sources. The Propagation Database is a key component of the operational tsunami forecasting system. At present, it includes a great number of pre-computed model runs for simulated earthquake sources in the Pacific, Atlantic and Indian Oceans (Titov, 2009; see also http://nctr.pmel.noaa.gov/propagation-database.html).

In the present study we suggest a simple method for short-term tsunami forecasting that does not require pre-computed tsunami waveforms. The method is based on a reciprocity principle (c.f. Rayleigh, 1945) and uses the only seismological information about co-ordinates of the earthquake epicentre and in situ tsunami data from one of the far-field open-ocean DART stations. This method can work in a real-time mode and, therefore, be effectively used for local tsunami warnings.

The method for a short-term tsunami forecast is based on Green's function for the wave equation possessing the fundamental property of symmetry. This property is well-known in acoustics and seismology as the reciprocity principle (Rayleigh, 1945). This principle is known as an asymptotic assumption.

The first use of this principle for long water waves was probably suggested by Loomis (Loomis, 1979). The reciprocity principle was found to be effective to optimize tsunami pre-computations. However, the effects of the ocean depth inhomogeneity had not been analysed previously. In view of large horizontal scales of tsunami sources, these effects can be significant. As a result of a study of the reciprocity principle and asymptotic conditions, the similarity conditions for reciprocal sources of long waves were obtained. Some applications of this principle to the tsunami research are considered in the current study. The numerical simulation was made using actual bathymetry. It was shown in numerical experiments that, with a quality sufficient for practical applications, the reciprocity principle works quite well, even when asymptotic conditions are not totally satisfied (i.e., when the characteristic time of the wave propagation between the reciprocity sources is only three times greater than the typical wave period and when the characteristic wave length is comparable with the size of the source area).

This paper presents a method of short-term tsunami forecasting using open-ocean sea level data from distant sites based on the reciprocity principle. The simple relationship (estimated transfer functions) enabled us to simulate tsunami waveforms for any selected oceanic point based only on the source location and sea level data from a remote reference site (Korolev, 2004, 2005). Information about the mechanism of the earthquake, the size of the centre or earthquake magnitude is not required.

This relationship (transfer function) consists of Laplace (Fourier) transforms. To obtain the target result, the inverse Laplace (Fourier) transform must be made. This problem is an ill-posed one. Two methods of inverse Laplace transform are described.

The method enables the transfer function to be evaluated during the event immediately after obtaining information about the earthquake epicentre location. The method enables us to compute expected tsunami waveforms in real-time mode and for any given ocean site. In contrast to the PMEL/NOAA technology, which is based on a pre-computed Propagation Database that needs updating for newly deployed open-ocean stations, the present method does not require pre-computed tsunami waveforms. The important advantage of this method is that it can be used irrespective of the actual source mechanism (seismic, submarine landslide or other phenomena).

The method was successfully applied to hindcast several recent tsunamis observed in the Pacific. The locations of the earthquake epicentres and the tsunami records from one of the NOAA DART sites were used as inputs for the modelling, while tsunami observations at other DART sites were used to verify and to estimate the effectiveness of this model.

The computed and observed tsunami waveforms for the regions of the Kuril Islands, Aleutian Islands and the US West Coast were in good agreement, satisfying the requirements of tsunami warning services (Korolev, 2011a, 2011b).

The good agreement is not only at points in open-ocean where DART stations are located, but also at points close to the coast of the Kuril Islands.

The proposed method can be considered as the basis for creating a program package that can be applied for early tsunami warnings for sites exposed to the threat of tsunamis. The method can be used for both regional and local tsunami warning services having access to the open-ocean (DART) data in a real-time mode. The result is a computational basis for decision making about issuing (confirmation) or cancelling a tsunami alarm. It is important that the method can be applied at newly developing tsunami warning centres.

2. Reciprocity principle

The reciprocity principle is a corollary of the properties of symmetry of Green's function for a linear wave equation. The Green's function is created, as the response of a medium on point-like (Dirac δ-function) disturbance. The reciprocity principle (relationship) expresses the same property of symmetry for sources with limited sizes. There are two aspects of this principle. A kinematic reciprocity principle says that the propagation time of a signal from point A to reciprocal point B is equal to the propagation time of a signal from point B to A and a dynamic reciprocity principle runs as follows: amplitudes of harmonic signals at points B (A) are identical, if the same source is located at A (B) (c.f. Rayleigh, 1945).

The validity of the reciprocity principle for an inhomogeneous medium with reflecting, absorbing and mixed (impedance) boundaries was proved later. It is valid for any area, including shadow zones and including the case of non-stationary waves (Chertock, 1970). The reciprocity principle is widely applied in hydroacoustics for hydrophone calibration (Uric, 1975).

Let us assume that a wave’s process is described by the linear wave equation. Let, in the space limited by reflecting and/or absorbing boundaries, the harmonic Source A be located at point A with co-ordinates ($$x_A,y_A$$), while the analogous harmonic Source B be at point B ($$x_B,y_B$$). Let also the pressure amplitude on the radiating surface be p_A, Q_A=p_AS_A, where S_A is the square of the radiating surface and p_A(B) is pressure amplitude in point B. The same is true for Source B: the pressure amplitude on the radiating surface is p_B, Q_B=p_BS_B, where S_B is the square of the radiating surface and p_B(A) is pressure amplitude at point A.

So the following reciprocity relationship (asymptotic equality) is valid (Brekhovskikh, 1960; Landau & Lifshitz, 1975)

$$p_A\left(B\right)·Q_B=p_B(A)·Q_A.$$E1

The obtained relationship is fair, if the following asymptotic conditions are valid:

1. the distance between the sources has to be significantly longer than the wave lengths;

2. horizontal scales of both sources areas have to be significantly shorter than the wave lengths.

In the case of an inhomogeneous medium, these conditions may be formulated in time terms:

a' the time of wave propagation between the sources should be significantly longer than typical wave periods;

b' the wave period should be significantly longer than typical time in the source.

Any connections of a source's parameters with parameters of a medium are not analysed. If for the acoustic processes, for example at ocean, the medium inhomogeneities are insignificant (variations of wave speed do not exceed 7%), then for tsunami waves they can give significant effects, as the speed of long waves differ by more than ten times (from 10 m/s on water depth 10 m up to 200 m/s on depth 4000 м). Considering large horizontal sizes of sources of a tsunami, the influence of these effects can be significant. Besides, the condition (b) can be defaulted, as the characteristic lengths of tsunami waves do not always exceed characteristic sizes of the centre. The problem arises as to whether or not the reciprocity principle for tsunami waves is valid. If yes, what are the conditions of validity and do they differ from conditions in acoustics? The conditions on sources should be investigated.

2.1. The reciprocity principle for non-stationary long waves

The deducing of the reciprocity relationship for tsunami waves is indicated in this paper with the purpose of showing analogies and differences from a traditional relationship. We follow the algorithm described in (Brekhovskikh, 1960; Landau & Lifshitz, 1975).

Let us assume that a wave’s process is described by the system of linear shallow-water equations, or wave equation, with variable water depth (Stoker, 1957)

$$\frac{{\partial }^2\zeta }{\partial t^2}-g\overrightarrow{\nabla }\cdot (D(x,y)\cdot \overrightarrow{\nabla }\zeta )=0$$E2

where ζ=ζ(t,x,y) – waveform, D(x,y) – water depth, g – gravity acceleration,$$\overrightarrow{\nabla }=(\frac{\partial }{\partial x},\frac{\partial }{\partial y})$$– differential operator.

Let Source A be located within area S_A with the centre at point A with co-ordinates x_A, y_A, while Source B be within area $$S_B$$ with the epicentre at B (co-ordinates x_B, y_B). For simplicity, the initial sea surface elevation is used as the wave source. ζ(t,x,y) is a waveform from Source A, η(t,x,y) is a waveform from Source B.

Initial conditions (t=0)

$$\zeta (0,x,y)={\zeta }_0(x,y)\text{in area}{S}_A$$E3

$$\eta (0,x,y)={\eta }_0(x,y)\text{in area}{S}_B$$E4

Boundary conditions

$$\zeta \left(t,x,y\right)=0\text{on}{\Gamma }_1\left(\text{absorbing boundary}\right)\text{,}$$E5

$$\zeta \left(t,x,y\right)=0 \text{on infinity}{\Gamma }_{\infty },$$E6

$$u_n\left(t,x,y\right)=0\text{on}{\Gamma }_2\left(\text{reflecting boundary}\right)\text{,}$$E7

$$\alpha \left(x,y\right)u_n\left(t,x,y\right)+\beta \left(x,y\right)\zeta \left(t,x,y\right)\text{on}{\Gamma }_3\left(\text{mixed, or impedance, condition}\right).$$E8

Here u_n(t,x,y) – normal component of mass velocity.

The same boundary conditions are applied to waveform η and mass velocity w_n from Source B.

While tsunamis are non-stationary waves, Laplace transform is applied to equation (2) and boundary conditions (5) for a waveform from Source A

$$s^2\zeta (s,x,y)-g\overrightarrow{\nabla }(D(x,y)\cdot \overrightarrow{\nabla }\zeta (s,x,y))=s{\zeta }_0(x,y).$$E9

The same is applied to waveform from Source B

$$s^2\eta (s,x,y)-g\overrightarrow{\nabla }(D(x,y)\cdot \overrightarrow{\nabla }\eta (s,x,y))=s{\eta }_0(x,y)$$E10

Multiplying the last but one equation by η(s,x,y) and the last equation by ζ(s,x,y) and subtracting one from another we obtain

$$\begin{array}{l}g\zeta (s,x,y)\cdot \overrightarrow{\nabla }(D(x,y)\cdot \overrightarrow{\nabla }\eta (s,x,y))-g\eta (s,x,y)\cdot \overrightarrow{\nabla }(D(x,y)\cdot \overrightarrow{\nabla }\zeta (s,x,y))=\\ =s(\eta (s,x,y)\cdot {\zeta }_0(x,y)-\zeta (s,x,y)\cdot {\eta }_0(x,y))\end{array}$$E11

This equation may be transformed to divergent form using relationship

$$f_1(x,y)\cdot \overrightarrow{\nabla }(D(x,y)\cdot \overrightarrow{\nabla }{f}_2(x,y))=\overrightarrow{\nabla }(D(x,y)\cdot f_1(x,y)\cdot \overrightarrow{\nabla }{f}_2(x,y))-D(x,y)\cdot \overrightarrow{\nabla }{f}_1(x,y)\cdot \overrightarrow{\nabla }{f}_2(x,y).$$E12

This yield

$$\begin{array}{l}g\cdot \overrightarrow{\nabla }(D(x,y)\cdot \zeta (s,x,y)\cdot \overrightarrow{\nabla }\eta (s,x,y))-g\cdot \overrightarrow{\nabla }(D(x,y)\cdot \eta (s,x,y)\cdot \overrightarrow{\nabla }\zeta (s,x,y))=\\ =s\cdot (\eta (s,x,y)\cdot {\zeta }_0(x,y)-\zeta (s,x,y)\cdot {\eta }_0(x,y))\end{array}$$E13

Then the integral is taken over space with boundaries Γ=Γ_∞+Γ₁+Γ₂+Γ₃. The integral of the left part of (6) is transformed into the integral along a contour Γ=Γ_∞+Γ₁+Γ₂+Γ₃ (Korn & Korn, 1961):

$$\begin{array}{l}g\cdot {\displaystyle \underset{\Gamma }{\oint }D\cdot \left[\zeta \cdot \left(\frac{\partial \eta }{\partial x}dy-\frac{\partial \eta }{\partial y}dx\right)-\eta \cdot \left(\frac{\partial \zeta }{\partial x}dy-\frac{\partial \zeta }{\partial y}dx\right)\right]}=\\ =-s\cdot {\displaystyle \underset{\Gamma }{\oint }D\cdot \left(\zeta \cdot w_n-\eta \cdot u_n\right)\cdot dl}\end{array}$$E14

The obtained integral over Γ_∞ and Г₁ is equal to 0 by virtue of the first and second conditions (5), the integral over Г₂ is equal to 0 by virtue of the third condition (5). The integrand is transformed to impedance form and the integral over Г₃ is equal to 0 by virtue of the fourth condition (5)

$$-s\cdot {\displaystyle \underset{{\Gamma }_3}{\int }D\cdot (\zeta \cdot w_n-\eta \cdot u_n)\cdot dl=}$$

$$-s\cdot {\displaystyle \underset{{\Gamma }_3}{\int }\frac{D}{\beta (x,y)}\cdot \left[\beta \cdot \zeta \cdot w_n+\alpha \cdot u_n\cdot w_n-\alpha \cdot w_n\cdot u_n-\beta \cdot \eta \cdot u_n\right]\cdot dl=}$$E15

$$=-s\cdot {\displaystyle \underset{{\Gamma }_3}{\int }\frac{D}{\beta }\cdot \left[\left(\alpha \cdot u_n+\beta \cdot \zeta \right)w_n-\left(\alpha \cdot w_n+\beta \cdot \eta \right)u_n\right]\cdot dl=0}$$E16

.

So, the integral of the left part of (6) is equal to 0. The integral of the right part (6) gives

$${\displaystyle \underset{S_A}{\int }\eta (s,x,y)\cdot {\zeta }_0(x,y)\cdot dS={\displaystyle \underset{S_B}{\int }\zeta (s,x,y)\cdot {\eta }_0(x,y)\cdot dS}}$$E17

Until now any suppositions, simplifying the problem, have not been done. Relationship (7) is the exact one.

The classical statement of the reciprocity principle requires the following conditions to be satisfied:

1. the distance between the sources has to be longer than the wave lengths (x, y scales of the source areas), i.e., the variation Δζ should be significantly smaller than function ζ(s,x,y) within the area S_B : Δζ<<ζ(s,x,y). This should also be true for the function η(s,x,y);

a' an alternative asymptotic condition may be formulated in the time domain: typical wave periods should be significantly smaller than the time of wave propagation between the sources;

1. the horizontal (x, y) scales of both source areas have to be shorter than the wave lengths of the generated tsunami waves.

If conditions (a), (a') and (b) are satisfied, functions ζ(s,x,y) and η(s,x,y) in (7) may be removed from the integrals and replaced by mean values at points B and A. Thus, an approximate (asymptotic) expression that follows from (7) is

$$\eta (s,x_A,y_A)\cdot {\displaystyle \underset{S_A}{\int }{\zeta }_0(x,y)\cdot dxdy}\approx \zeta (s,x_B,y_B)\cdot {\displaystyle \underset{S_T}{\int }{\eta }_0(x,y)\cdot dxdy}$$E18

.

Designating $${\displaystyle \underset{S_A}{\int }{\zeta }_0(x,y)\cdot dxdy}=Q_A$$ and$${\displaystyle \underset{S_B}{\int }{\eta }_0(x,y)\cdot dxdy}=Q_B$$, where $$Q_A$$ and $$Q_B$$ are the volumes of the initial disturbances, we can rewrite the last relationship as

$$\eta (s,x_A,y_A)\cdot Q_A\approx \zeta (s,x_B,y_B)\cdot Q_B$$E19

The form of relationship (8) matches with the form of traditional relation (1).

Due to the medium inhomogeneity (variable speed of wave propagation) and large horizontal source scales in comparison with the water depth, the initial free surface elevations should satisfy the following conditions based on common reasons of similarity and dimensions:

• the forms of the initial disturbances should be similar;

• the disturbance radii, $$R_A$$and$$R_B$$, are related as

$$\frac{R_A}{\sqrt{g{D}_A}}=\frac{R_B}{\sqrt{g{D}_B}},$$E20

where D_A and D_B are the water depths at points A and B, and g is the gravitational acceleration.

The solution of one problem confirms this relationship (see Appendix A).

The approximate (asymptotic) relationship (8) and the condition of similarity (9) is the dynamic reciprocity principle (relationship) for non-stationary long water waves (Poplavsky et al., 1997).

Additionally, if the disturbance volumes are identical (a_AR_A²=a_TR_T², i.e. Q_A=Q_B), then after the inverse Laplace transform, the waveforms will also match each other

$${\zeta }_T(t,x_A,y_A)={\zeta }_A(t,x_T,y_T).$$E21

2.2. The analysis of the reciprocity principle for long waves

It is of interest to identify whether relationship (10) is fair if conditions (a) and (b) are not satisfied, therefore, this will be the subject of the following study.

To verify this approach we made several numerical experiments. The numerical simulations based on actual bathymetry were conducted for the region of the Southern Kuril Islands. Source 1 was located in the Sea of Okhotsk northwest of Kunashir Island; Source 2 was located in the Yuzhno-Kurilsk Strait, between the Kunashir and Shikotan Islands; and Source T was located in the Pacific Ocean southeast of Shikotan Island. The computational domain for the numerical experiments indicating the location of the sources is shown in Fig. 1. The bathymetry grid for these experiments was chosen to have a spatial step $$\Delta x=\Delta y$$= 1000 m.

The parameters of the initial free surface elevations were chosen taking into account the above-mentioned conditions. They were the following:

Source 1: Max ampl. a₁ = 0.4 m, radius R₁ = 25 km, water depth in epicentre D₁= 2418 m;

Source 2: Max ampl. a₂ = 15.6 m, radius R₂ = 4 km, water depth in epicentre D₂ 63.6 m;

Source T: Max ampl. a_T = 0.49 m, radius R_T = 22.5 km, water depth in epicentre D_T= 2030 m.

The results of these simulations are shown in Fig. 2. The correlation coefficient, ρ, was chosen as a criterion for comparison of the reciprocal waveforms.

Figure 1.

A map showing the region of the South Kuril Islands and a scheme of numerical experiments. Symbols 1, 2 and T indicate the epicentres of the sources. Initial sea surface elevations at Source 1 and Source T are pictured.

The coincidence of waveforms according to (10) is quite good despite wave sheltering by a chain of islands. The typical wave periods are about 10-18 min; the time of wave propagation is 57 min between Sources 1 and T, 50 min between Sources 1 and 2, and 40 min between Sources 2 and T. It is possible to make certain that wave lengths in the reciprocity source areas slightly exceed their sizes. For example, the wave length of the wave from Source 1 (λ₁₂) at Site 2 for the period of T = 18 min and water depth D₂ = 63.6 m is $${\lambda }_{12}=T\cdot \sqrt{g{D}_2}\approx$$ 27 km, while the diameter of Source 2 is 8 km.

Figure 2.

Illustration of the reciprocity relationships. (a) Waveforms at Site T due to Source 1 (blue solid line) and at Site 1 due to Source T (red dashed line); correlation coefficient ρ_1T = 0.96. (b) Waveforms at Site 2 due to Source 1 (blue solid line) and at Site 1 due to Source 2 (red dashed line); correlation coefficient ρ₁₂ = 0.94. (c) Waveforms at Site T due to Source 2 (blue solid) and at Site 2 due to Source T (red dashed line) correlation coefficient ρ_2T = 0.91.

The wave length of the wave from Source T (λ_T2) at Site 2 for period of T = 8 min and water depth D₂ = 63.6 m is $${\lambda }_{T2}=T\cdot \sqrt{g{D}_2}\approx$$ 12 km, while the diameter of Source 2 is 8 km.

Numerical experiments show that relationship (8) with condition (9) is true even if the time of wave propagation between sources is only 3 – 4 times longer than wave periods and if the wave lengths are comparable to or slightly longer than these sizes. As was found from our experiments, changing the radius of one of the sources by 25% does not significantly influence the output.

Additional experiments were done when Sources 1 and 2 were the same, but Source T was not a circle but an ellipse. The ellipse minor ($$b_2$$) and major ($$b_1$$) axes were chosen as (first) $$b_1$$= 63.6 km and $$b_2$$= 31.8 km ($$b_2/b_1$$= 1/2); and (second) $$b_1$$= 77.9 km and $$b_2$$= 26.0 km ($$b_2/b_1$$= 1/3); maximum amplitude a_T = 0.49 m and water depth at epicentre D_T = 2030 m (for both cases). We also changed the orientation of the major axis relative to true north. The condition of similarity of the source forms was not true, but it was accepted that $$\frac{2R_{1(2)}}{\sqrt{g{D}_{1(2)}}}=\frac{\sqrt{b_1\cdot b_2}}{\sqrt{g{D}_T}}$$ and$$4a_1{R}_1^2=a_T{b}_1{b}_2$$. The respective numerical results are shown in Fig. 3.

Figure 3.

Fig. 3. Illustration of the reciprocity relationships. Elliptic Source T is elongated to the northwest (see Fig.1). The ratio of minor/major axes is equal to 1/2 (left column) and to 1/3 (right column). Upper row: waveforms (a) at Site T due to Source 1 (blue solid line), correlation coefficients are ρ_1T = 0.90 and (b) at Site 1 due to Source T (red dashed line), ρ_T1 = 0.87. Lower row: waveforms (c) at Site T due to Source 2 (blue solid line), correlation coefficients are ρ_2T = 0.86 and (d) at Site 2 due to Source T (red dashed line), ρ_T2 = 0.78.

The correlation coefficients in all experiments are in the range of 0.67 to 0.90. The best result is in the case when the major axis was directed toward the northwest. This agreement is good enough and sufficient for any practical applications. Relationship (10) was found to be fair not only for circular sources, but also for elliptic sources, at least for the moderate ratios of axes' lengths.

3. The method of short-term tsunami forecasting

The classical statement of the reciprocity principle requires conditions (a) and (b) (see the previous section), however, our numerical experiments demonstrate that it can still be applied without these strict contingencies. The results presented below could become the basis for the short-term tsunami forecasting.

3.1. Derivation of relationships for a short-term tsunami forecast

It was demonstrated (Korolev, 2004, 2005) that the derived reciprocity relationships could be efficiently applied for the problem of short-term (operative) tsunami forecasts based on sea level data from a remote open-ocean site. It is assumed that reciprocity relationship (8) with similarity condition (9) is true when one of the sources is natural, while the other source is artificial (numerical).

Let's tsunami Source T, the initial sea surface elevation, is located within the area S_T. The total volume of the disturbed uplift is Q_T. Let us also assume that M is the site of a sea level gauge (i.e., DART site) and A is the site of interest (a point for tsunami forecasting). The problem is to compute the tsunami waveform at site A (“target site”) using the sea level data at site M (“reference site”). Mathematically, this problem does not have a unique solution, however, it can be considered as a "consumer problem”. Specific methods should be developed for solving this problem.

In this section, “waveform” means the Laplace (Fourier) transform (spectra) of the waveform. Let the reciprocal sources be located within areas with epicentres at A and M. Their disturbed volumes are Q_A and Q_M, while η_A(s,x_T,y_T) and η_M(s,x_T,y_T) are the waveforms at site T from reciprocal sources S_A and S_M. For the natural tsunami waveform ζ_T from tsunami Source T and for waveforms from the reciprocal sources in S_A and S_M, relationships (8) are

$$\begin{array}{l}{\zeta }_T(s,A)\cdot Q_A={\eta }_A(s,T)\cdot Q_T;\\ {\zeta }_T(s,M)\cdot Q_M={\eta }_M(s,T)\cdot Q_T.\end{array}$$E22

.

In these expressions we designated: ζ_T(s,A)≡ζ_T(s,x_A,y_A), ζ_T(s,M)≡ζ_T(s,x_M,y_M) and so on.

In addition, let's assume that the auxiliary sources are in the areas with the same epicentres at A, M and T (waveforms and disturbed volumes are designated by symbols with primes). For waveforms from these sources the reciprocity relationships are also true

$$\begin{array}{l}\eta {\text{'}}_T(s,A)\cdot Q{\text{'}}_A=\eta {\text{'}}_A(s,T)\cdot Q{\text{'}}_T\\ \eta {\text{'}}_T(s,M)\cdot Q{\text{'}}_M=\eta {\text{'}}_M(s,T)\cdot Q{\text{'}}_T\end{array}$$E23

As the auxiliary sources may be selected rather arbitrarily, in particular, Q_A, Q_M and Q'_A, Q'_M may be selected in such a way that η'_A is equal to η_A and similarly for η'_M and η_M. It is true if Q_A = Q'_A and the same for Source M.

Eliminating unknown Q_T from the above-mentioned two sets of equations yields the following relationship

$$\frac{{\zeta }_T(s,A)}{{\zeta }_T(s,M)}=\frac{{\eta }_T(s,A)}{{\eta }_T(s,M)}$$E24

Let us assume that functions on the left side of the equation relate to one tsunami source and functions on the right side relate to another source, and that both sources have the same epicentre. Relationship (11) supports the well-known fact that spectra of any two tsunamis from different sources are similar when recorded at the same site, whereas spectra of any two tsunamis from the same source are dissimilar at different sites (Takahasi & Aida, 1961; Miller, 1972).

From (11) we can get the following main relationship

$${\zeta }_T(s,A)={\zeta }_T(s,M)\frac{{\eta }_T(s,A)}{{\eta }_T(s,M)}$$E25

Here ζ_T(s,A) is the target function (tsunami waveform at A); ζ_T(s,M) is the reference function (tsunami waveform at M, i.e., open-ocean DART data); and η_T(s.M) and η_T(s,A) are the numerical waveforms at M and A from the auxiliary source located at S_T. All functions are Laplace transforms. The ratio on the right side of (12) plays the role of a transfer function. The transfer function allows computing a tsunami waveform at a target point on ocean level data at a reference site.

As the sea level data and computed data are digital ones, so the functions included in (12) are represented as numerical Laplace transformations or z-transformation, for example $$\zeta (s,A)=\zeta (z,A)={\displaystyle \sum _{k=0}^N{\zeta }_k(A)e^{-ks\tau }}={\displaystyle \sum _{k=0}^N{\zeta }_k(A)z^k}$$ (see Appendix B), where ζ_k(A) – is the value of function ζ(t,A) for a moment t=kτ, τ – is a time sampling. Due to this, (12) may be rewritten as

$$\zeta (z,A)={\displaystyle \sum _{k=0}^N{\zeta }_k(A)\cdot z^k}={\displaystyle \sum _{k=0}^{N_1}{\zeta }_k(M)\cdot z^k}\frac{{\displaystyle \sum _{k=0}^{N_2}{\eta }_k(A)\cdot z^k}}{{\displaystyle \sum _{k=0}^{N_3}{\eta }_k(M)\cdot z^k}}$$E26

The arrival time of the forecasted tsunami T_A at target point A can be estimated as T_A = T₁ + T₂-T₃, where T₁ - time of tsunami propagation to the point of registration M, T₂ - an arrival time of wave from an auxiliary source to point A, T₃ - an arrival time of wave from an auxiliary source to point M.

At the first stage, while the record obtained from the station DART has a short duration (for example, one half period of the first wave, N₁ = 10-20 minutes), the lengths of all power series in (13) should be equal: N₁ = N₂ = N₃ = N, otherwise the result will be inadequate. But the approximate forecast is possible if N₃ = N₁, but N₂ is significantly more than N₁. Such a forecast is approximate, but probably takes into account secondary waves. As the sea level information from the DART station at M is received, the duration and quality of the forecast is increased. The computations with variable duration of level data allow making the forecast in real-time mode (see Appendix B).

The practical interest is to predict tsunami waveforms of long duration adequately, which would take into account secondary waves reflected from island chains, submarine ridges and mainland coasts. These secondary waves have amplitudes which often significantly exceed amplitudes of first tsunami waves.

Such a tsunami forecast of long duration can be made on the limited information about a level. Often the tsunami in the ocean only has some first waves with significant amplitudes. If we neglect waves with small amplitude in a wave tail, the duration N₁ should be equal to the duration of the registered tsunami with significant waves and the duration N₂ could somehow be long, but durations N₃ and N₁ should be equal. The forecast duration N in this case is equal to N₂. In this way secondary waves can be taken into account.

The moment of the preliminary forecast completion is determined by the time of tsunami propagation up to a point of registration and the time necessary for tsunami identification (the first half period of the first wave).

The division and the multiplication of polynomials in the right part of (13) could give the required result. Unfortunately, the direct application of the offered algorithm, gives, as a rule, a diverging result. This problem is an ill-posed one. Two methods of the solving this problem (regularization methods) permit to obtain a converging result and are described in Appendix B.

It should be noted that the generation mechanism of the actual tsunami is not important here. The method can be applied for short-term forecasting irrespective of the tsunami source, in particular for tsunamis produced by earthquakes, volcanic eruptions, submarine landslides or other phenomena. However, for the practical realization of the proposed method it is important to assume that the transfer function for the corresponding event can be constructed as the ratio of the Laplace transforms of the waveforms at appropriate sites generated by an initial free surface elevation. This elevation was supposed to be circular with a diameter of 50-100 km which is the characteristic transversal size of the earthquake source. There is no requirement to know the detailed seismological information about the source; it is enough to know the co-ordinates of the earthquake (or another source) epicentre. This will enable us, with accuracy sufficient for practical use, to create a short-term tsunami forecast for specific coastal sites.

In fact, the proposed method is a generalization of the methods suggested earlier ( Chubarov & Shokin, 1995; Chung et al., 1995).

3.2. Dependence of the hindcast quality on the diameter of the auxiliary source

There are questions about how the size of auxiliary sources affects the results of the computations.

Numerical experiments with auxiliary sources of different diameters (20, 30, 40, 50, 75 and 100 km) for sites in the ocean gave almost identical results.

3.3. Computing the short-term tsunami forecast

The algorithm for computing the tsunami waveform based on sea level data from remote sites includes the following steps:

1. deriving the co-ordinates of the earthquake (or other source) epicentre;

2. computing the tsunami waveforms for the reference site M and the target (forecasting) site A, i.e., functions η_T(s,M) and η_T(s,A). The auxiliary source is a circle with diameter of 50-100 km located at the centre of the actual tsunami source. The computation (evaluation of the transfer function) should be completed before the tsunami reaches the reference DART site (normally the site nearest to the source);

3. acquiring to the tsunami data from the appropriate reference DART station; eliminating tidal components from the corresponding record;

4. after recording and identifying tsunami waves at the reference DART station (Site M), the tsunami forecast can be made for the target site A based on relationship (12, 13). Similarly, in real-time mode, tsunami forecasts can be provided for other target sites (e.g., for particular sites along the coast). In fact, the preliminary forecast can be given after recording the first tsunami semi-wave. Hereafter, as new and more complete information is received, this forecast can be corrected and improved;

5. the forecast results can be used by TWSs to assist in making their decision about declaring or cancelling the tsunami warning (alarm).

4. Application of the method for short-term tsunami forecasting

Some preliminary research results published in 2005 (Korolev, 2005) demonstrated that the proposed method can be efficiently used for short-term tsunami forecasting. In particular, it was shown that the proposed method can work in a real-time mode and that the errors in the epicentre location only slightly affect the quality of the forecast.

Some new results were published recently by Korolev (2011a, 2011b) and some of these are represented below.

4.1. Preliminary information

In the present study we apply the method of short-term tsunami forecasting to hindcast tsunami waveforms for some tsunamigenic earthquakes. Earthquakes occurred in 2006, 2007 and 2009 in the region of the Central Kuril Islands, eastward of Simushir Island, near the Chilean coast in 2010 and eastward of Honshu Island in 2011. Associated tsunamis were recorded by some tide gauges on the coasts of Russia and by a number of DART stations located along the Kuril Islands, Japan, the Aleutian Islands and the US West Coast (NDBC DART, n.d.).

Tsunami waveforms were simulated at the locations of the DART stations and near to the coast of Kuril Islands. The actual in situ tsunami data from the DART stations closest to the earthquake epicentres were used for the hindcast. The data from other DART stations, located further from the epicentre, were used to compare with the predictions and to verify the method. The correlation coefficient was chosen as criterion of coincidence. In the present study we used deep-ocean stations not affected by coastal resonant effects and surf beats. All computations were made based on the shallow-water numerical model (Poplavsky et al., 1997). We used the bathymetric data that were interpolated from the ETOPO2 global dataset (Smith & Sandwell, 1994).

No seismological information about the earthquake, except the epicentre co-ordinates, was used in the following computations.

The scheme of the research experiments, indicating the earthquake epicentre and positions of the DART stations, is shown in Fig. 4. The hindcasting results are represented below.

4.2. Simushir tsunamis of 2006, 2007 and 2009

The auxiliary source with an initial circular sea surface uplift of 75 km in diameter and with a maximum height of 10 m was chosen for transfer function creation in the cases of the 2006 and 2007 events. The selected source centre coincided with the earthquake epicentre.

Results were published in (Korolev, 2011b).

Figure 4.

A map of the Pacific Ocean showing positions of DART stations (from NDBC DART, n.d.) and epicentres of some earthquakes.Digits denote: 1 – 2006 Simushir tsunami; 2 – 2007 Simushir tsunami; 3 – 2009 Simushir tsunami; 4 – 2010 Chile tsunami; 5 – 09.03.2011 tsunami; 6 – 2011 Great Tohoku tsunami.

4.2.1. The November 2006 Simushir tsunami

The tsunami was generated by a major earthquake (M_w = 8.3) on November 15, 2006 at 11:14 UTC with the epicentre at 46.592ºN, 153.266ºE near the Central Kuril Islands (NOAA/WDC, n.d). Significant wave heights were observed along the coasts of the Kuril Islands (c.f. Rabinovich et al., 2008; McInnes et al., 2009).

The November 15, 2006 Simushir tsunami was accompanied by a big run-up at islands nearest to the earthquake centre, but small amplitudes in the northern and southern Kuril Islands, and at the same time caused considerable damage in Crescent City on the US West Coast (Horrillo et al., 2008; Kowalik et al., 2008).

A tsunami alarm was declared at 11:30 UTC by Sakhalin Tsunami Warning Center on the coast of all Kuril Islands and an evacuation of people from dangerous areas was done. Vessels were ordered to sail out into open sea. The alarm was cancelled at 13:19 UTC.

The proposed method was applied to simulate the operational tsunami forecast in the North Pacific Ocean along the Aleutian Islands and the US West Coast. To compute the tsunami waveforms, the data of DART stations 21414 and 46408 were used as the basis (see Fig. 4 for the DART positions). The ocean depth at the epicentre, based on the bathymetry grid readings, was 2803 m. Spherical co-ordinates with a spatial grid step of 3.8 km at latitude 40ºN were used.

Figure 5 shows hindcast results and the waveforms of the observed tsunami waves.

Figure 5.

Computed and observed tsunami waveforms for the 2006 Simushir earthquake for DART sites in the North Pacific Ocean. Left column: using DART 21414 as the reference station; right column: using DART 46408 as the reference station.

The computed waveforms match well the observed tsunami waveforms at DART stations 46402 and 46409 along the Aleutian Islands. Similar good agreement was observed also for DART stations 46419 and 46411 along the US West Coast. In all cases, the best agreement was for a few leading waves; the correlation coefficients for these waves were within the range $$\rho$$ = 0.70–0.85.

4.2.2. The January 2007 Simushir tsunami

The tsunami was generated by a major earthquake (M_w = 8.1) on January 13, 2007 at 4:23 UTC with the epicentre at 46.243ºN, 154.524ºE located on the oceanic slope of the Central Kuril Trench (NOAA/WDC, n.d.). This tsunami was weaker than the 2006 Simushir tsunami, but was locally quite intensive (Rabinovich et al., 2008).

The tsunami alarm was declared on the coast of the Kuril Islands at 4:36 UTC.

Tsunami waves were recorded with amplitudes at Yuzhno-Kurilsk 5 - 6 cm and at Malokurilskoe 36 cm. Due to the weak manifestation of the tsunami on the coasts of the Kuril Islands, Kamchatka and the Aleutian Islands, the tsunami alarm was cancelled at 7:49 UTC.

This tsunami was recorded by a smaller number of DART stations than the 2006 tsunami (c.f. Laverov et al., 2009), however, one of these stations was DART 21413, which was not in operation during the previous (2006) event. This station was located southward from the 2007 epicentre (690 nautical miles southeast from Tokyo). During the 2006 event the tsunami was only recorded by DART stations located to the east of the source. It is of additional interest to compare the results for open-ocean stations located in various directions from the source. Simulation of the operational tsunami forecast was done. The data of DART stations 21414 (east of the source) and 21413 (south of the source) were used as the basis (see Fig. 4 for the DART positions). The ocean depth at the epicentre was 6877 m. Spherical co-ordinates with a spatial grid step of 3.8 km at latitude 40ºN were used. The corresponding results are shown in Fig. 6.

In general, the agreement between computed and observed waveforms is satisfactory both for DART 46408 near the Aleutian Islands and for DARTs 46419 and 46412 along the US West Coast. Hindcasts based on data from both stations 21414 (eastward from the source) and 21413 (southward from the source) are of almost identical quality. For all cases, the computations show correctly the tsunami manifestation during the initial (negative) phase. The correlation coefficients were within the range $$\rho$$ = 0.50-0.72.

Figure 6.

The same as in Fig. 5, but for the 2007 Simushir tsunami. Left column is for DART 21414 as the reference station; right column is for DART 21413 as the reference station.

4.2.3. The January 2009 Simushir tsunami

This earthquake (M_w = 7.4) occurred on January 15, 2009 at 17:49 UTC approximately in the same region of the Central Kuril Islands as the 2006 and 2007 earthquakes, but it was weaker than the two previous earthquakes. The epicentre of the earthquake was located at 46.857ºN, 155.154ºE (NOAA/WDC, n.d.). The earthquake generated a tsunami that was recorded on the coast of the Kuril Islands and by a few DART stations, in particular by DART 21416 located 240 miles from the Kamchatka Peninsula. The latter station was the nearest to the epicentre. A tsunami alarm was not announced during this event.

Figure 7 shows a more detailed map than in Fig. 4, indicating epicentres of earthquakes, DART positions in the northwest Pacific and target points on the coast of the Kuril Islands.

Figure 7.

Map fragment of northwest Pacific. Conventional signs: S-K – Severo-Kurilsk, K - Kurilsk, B – Burevestnik (Kasatka Bay), Yu-K – Yuzhno-Kurilsk; asterisks – locations of earthquake epicentres (and dates), rhombuses – locations of DART stations.

The location of DART station 21416 is out of the scheme, the true position is to the east of that marked on the map at longitude 163° 30'E, the true site of DART station 21413 is to the south of that marked on the map at a point with latitude 30° 31'N.

The auxiliary source was taken as a circular uplift of the sea surface 50 km in diameter and maximum height of 8 m; the source centre coincided with the earthquake epicentre. This source was used as the input for the model. The ocean depth at the epicentre was 6645 m.

The computed and observed tsunami waveforms are shown on Fig. 8. The data from DART 21416 were used to hindcast the waveforms for DART 46408 and 46413, located eastward from the source, and for DART 21413, located southward from the source. Similar computations were also made for DARTs 46408, 46413 and 21416 based on the data from DART 21413.

Figure 8.

The same as in Fig. 5, but for the 2009 Simushir tsunami. Left column is for DART 21416 as the reference station; right column is for DART 21413 as the reference station.

The computed waveforms are in good agreement with the observed tsunami records for both the area of the Kamchatka Peninsula and the area of the Aleutian Islands, especially for the leading waves. The correlation coefficients were within the range $$\rho$$ = 0.75–0.80.

4.2.4. The 2009 Simushir tsunami on the Kuril Islands

Previous results demonstrated the possibilities of the proposed method. More interesting is the answers the following questions. How can the proposed method produce a tsunami forecast near to the coast? Can the proposed method predict secondary waves of large amplitudes? How well can the proposed method estimate the duration of the alarm period?

The auxiliary source to create a transfer function was an axisymmetric free surface elevation with a diameter of 50 km and a maximum amplitude of 8 m. Water depth at the earthquake centre based on the bathymetry grid readings was equal to 6887 m. Mercator projection with a spatial grid step 1.1 km was used with latitude of scale 52º.

The data of the closest to the earthquake centre DART station 21416, located north-east of the source, were used for the simulation of the tsunami forecast on the coast of the Kuril Islands. The sea level recorded by this station is shown in Fig. 9.

Figure 9.

The sea level recorded by DART station 21416.

The January 15 tsunami was recorded by a tide gauge at Severo-Kurilsk, located north of the source. The tsunami waveform was computed at the point near Severo-Kurilsk in the strait at a distance of about 1.5 km offshore and at a depth of 20 m.

Results of simulation are represented in Fig. 10. This figure shows a tide gauge record at Severo-Kurilsk with a duration of 17 hours. The specific natural oscillations of the strait up to 19 hours on the sea level record on January 15 can be seen. Then the nature of the oscillations changed due to the tsunami. After 6 - 7 hours the specific mode of natural oscillations was restored, but with increased amplitude.

The blue line designates the recorded form of oscillations of a sea level at the port of Severo-Kurilsk, red - hindcast:

1. preliminary forecast at the point S-K (see Fig. 8) based on data of DART station 21416 with a duration of 20 min;

2. updated forecast at the same point based on data of DART station 21416 with a duration of 82 min.

On Fig. 10 the vertical solid line shows the time of the main shock (17:49 UTC), the vertical dashed line shows the moment of the forecast.

The preliminary estimated tsunami waveform based on 20 minute data of DART station 21416 is shown in Fig. 10 (a). The updated forecast using 82 minute level data of DART station21416 is shown in Fig. 10 (b).

Forecast, preliminary and adjusted, shows that the tsunami arrival is expected at 19:20, 90 min after the main shock. Waves with a greater amplitude followed the first wave are expected 1.5 hours later. The duration of the tsunami is estimated to be 6 - 7 hours.

Figure 10.

Results of the 2009 Simushir tsunami hindcast close to Severo-Kurilsk.

The computations necessary to create a transfer function are performed for 15 min, while the tsunami travel time to the point 21416 equals to 41 min, to Severo-Kurilsk - 75 min. Taking into account the time of tsunami identification (20 min), the preliminary tsunami forecast at Severo-Kurilsk may be given 30 minutes before the first wave arrives at this point.

The above example shows the principal possibility of tsunami forecast near to the coast of the Kuril Islands based on data of level stations, located on the ocean side of the Kuril-Kamchatka Trench.

4.3. 2010 Chile tsunami

The tsunami was generated by earthquake (M_w = 8.8) which occurred on February 27, 2010 at 6:34 UTC near the Chilean coast. The co-ordinates of the earthquake epicentre are 36.122ºS, 72.898ºW (NOAA/WDC, n.d.). The tsunami was recorded by a number of DART stations in the Pacific.

4.3.1. Hindcast of 2010 Chile tsunami in ocean

Post-event computations simulating the short-term tsunami forecast are based on data closest to the epicentre of the earthquake recorded at DART station 32412 over 30 minute duration (Fig. 11). An auxiliary source to create the transfer function was defined as a circular initial free surface elevation with a diameter of 100 km and a maximum height of 10 m with the centre coincident with the epicentre of the earthquake. The computations were performed in spherical co-ordinates with a spatial grid step 5 km at the equator.

The results of hindcasting the tsunami waveforms in the Pacific (see Fig. 4 for the positions of DART stations and tsunami epicentre) are represented in Fig. 12.

The results of the simulation show a good agreement between hindcasted and observed tsunami waveforms in the north and the west directions from the earthquake centre.

Figure 11.

Sea level recorded by DART station 32412.

Figure 12.

Results of hindcasting the 2010 Chile tsunami in the Pacific based on DART 32412 data.

4.3.2. Hindcast of 2010 Chile tsunami near the coast of the Kuril Islands

The earthquake occurred on February 27 at 06:34 UTC. The tsunami alarm was declared by the Sakhalin Tsunami Warning Center on all the Kurile Islands on February 28 at 00:50 UTC (10:50 local time) 3 h 22 min before the first wave arrived at Severo-Kurilsk (Paramushir Island). The population of all settlements of the Kuril Islands was evacuated and ships were withdrawn to the open sea. The wave amplitudes for 4 hours did not exceed 0.3 m. The duration of the alarm was about 4 hours. Shortly after the alarm was cancelled (4 hours 8 minutes after the first wave arrival) a wave with an amplitude of 1 m was recorded. Once again the tsunami alarm had not been announced.

For tsunami prediction near the coast of the Kuril Islands, computations have been made from the source of the tsunami off the coast of Chile, initially using spherical co-ordinates with a spatial grid step 5 km at the equator. Then computations were made using a more detailed grid (left bottom corner 37° 51'N, 139° 26'E, right top corner 45° N, 180° E). The spatial grid step was 1 km on scale latitude 45°N. Part of the scheme is shown in Fig. 7.

The tsunami hindcasting is based on data of DART station 32412 with a 30 minute duration The hindcast result is represented in Fig. 13.

Figure 13.

Chile tsunami hindcast at Severo-Kurilsk (point S-K on Fig. 7) using data of DART station 32412.

The simulated tsunami waveform coincides well with the form of the observed tsunami. The wave amplitudes during the first 4 hours do not exceed 0.5 m. After 4 hours the appearance of waves with increased amplitude up to 1 m is correctly predicted.

4.4. Tohoku foreshock 2011

The weak tsunami M_w=7.2 which arose off east coast of Honshu Island on March 09, 2011 (epicentre co-ordinates 38.510°N, 142.792°E (NOAA/WDC, n.d.)) was recorded by a small number of DART stations. The tsunami hindcast using data of the 21418 DART station, nearest to the centre, was made in the area eastward of the Kuril Islands at points 21401, 21419 and in the Pacific Ocean at a point 21413 1300 km to the southeast of Tokyo (see Fig.7 for DART positions).

The auxiliary source was taken as a circular uplift of the sea surface 100 km in diameter and maximum height of 10 m; the source centre coincided with the earthquake epicentre. The ocean depth at the epicentre was 1374 m. We used the Mercator projection with a spatial grid step 5 km at latitude 45°N.

Represented in Fig. 14, the result shows a good coincidence between computed and actual tsunami waveforms.

The computations were executed for points in the open-ocean. The various edge, resonant effects do not influence tsunamis in the open-ocean. This result demonstrates the possibilities of the proposed method. Based on data of DART station 21418, located to the east of the earthquake epicentre, the tsunami waveforms are obtained to the northeast of the epicentre (points 21401 and 21419) and to the southeast (point 21413).

Figure 14.

Hindcast of the March 9, 2011 tsunami in the Pacific.

4.5. 2011 Great Tohoku tsunami

4.5.1. 2011 Great Tohoku earthquake and tsunami

This earthquake with magnitude M_w = 9.0 off the northeast coast of Honshu Island occurred on March 11 at 5:46 UTC. The co-ordinates of the epicentre are 38.322ºN, 142.369ºE (NOAA/WDC, n.d.).

As a result of the earthquake a tsunami occurred. It was anomalously high and very destructive on the coast of Honshu Island. Against the backdrop of a long rise in sea level with a height of 2 metres there was an additional shorter rise, with increased amplitude up to 5 m 11 minutes after (Fig. 15).

The great tsunami resulted from a combination of crustal deformations of the ocean floor due to up-thrust tectonic motions, augmented by additional uplift due to the quake's slow and long rupturing process, as well as large coseismic lateral movements which compressed and deformed the compacted sediments. This mechanism is discussed in (Pararas-Carayannis, 2011).

Figure 15.

Registration of the 11.03.2011 tsunami near to the centre. Positions of earthquake epicentre with М = 9.0 and sea level stations ТМ1 and ТМ2 - at the left. The form of tsunami registered by stations - on the right (see http://outreach.eri.u-tokyo.ac.jp/wordpress/wp-content/uploads/2011/03/TsunamiKamaishiMeter.jpg).

4.5.2. The hindcast of 2011 Great Tohoku tsunami in the ocean

A simulation of operative tsunami forecasting for March 11, 2011 was made for points in the open-ocean. The data of DART station 21401 were used in the computations. The relative level of the ocean registered by station 21401 is shown in Fig. 16.

Figure 16.

Sea level recorded by DART station 21401 on March 11, 2011 (the tidal components are deleted), the vertical red line – the moment of the main shock.

To create the transfer function the auxiliary computations of the waveform at the corresponding points from the source in the form of an initial circular elevation of the free-surface with diameter 100 km and an amplitude 10 m are made. The centre of the source coincides with the epicentre of the earthquake, the depth of the ocean according to the bathymetric grid equals to 887 m. Computations were made using spherical co-ordinates with the spatial grid step 5 km at a scale latitude 0°. The computations were based on a 30 minute tsunami duration of the data from 21401 DART station.

Figure 17 demonstrates the high quality of the coincidence of simulated and observed tsunami waveforms both in the northeastern and southern directions from the source. Durations exceeding 30 min in obtained data from DART station 21401 used in the computations do not affect the quality of the forecast.

Figure 17.

Hindcast of the March 11, 2011 Great Tohoku tsunami in the Pacific.

Analogous simulations have been performed using other means by Titov, Fujii and Satake (Japan (East Coast of Honshu) Tsunami, 2011; Fujii & Satake, 2011). The quality of the coincidence of simulated and observed tsunami waveforms in these works as well as in the presented work is comparable.

4.5.3. The hindcast of 2011 Great Tohoku tsunami near the coast of the Kuril Islands

This strong earthquake off the northeast coast of Honshu Island occurred on March 11 at 5:46 UTC.

The tsunami alarm was declared by the Sakhalin Tsunami Warning Center at 5:58 UTC (local time 15:58) on all Kuril Islands. The population deemed to be in the threatened areas was evacuated to safe places and vessels were recommended to sail out into the open sea.

At 6:43 deep-sea DART station 21401 (located southeast of Iturup Island) registered a tsunami wave with amplitude 67 cm.

The oscillations of sea level, caused by the tsunami, with amplitude of about 1 metre lasted until March 12, therefore, the duration of the tsunami alarm was about 20 hours. The tsunami alarm on all Kuril Islands was cancelled on March 12 at 2:20 UTC (local time 12:20).

The tsunami arrived at Severo-Kurilsk at 9:00. From information from vessels in the sea near Severo-Kurilsk the sea level under keel changed from 4.2 m up to 2.6 m, the swing was 1.6 m. Recorded wave amplitudes were about 0.3 m in Kurilsk. According to visual observations the amplitudes of the tsunami in Burevestnik were approximately 1 m. The tsunami reached Yuzhno-Kurilsk at 7:44, the recorded wave amplitudes were about 1 m. The tsunami alarm was declared 1.5 to 3 hours before the tsunami arrival at different points of the Kuril Islands.

To create the transfer function the auxiliary computations of the waveform at the corresponding points from the source in the form of an initial circular elevation of the free-surface with diameter 100 km and an amplitude 10 m were made. The sea depth at the epicentre is equal to 970.1 m. The centre of the source coincides with the epicentre of the earthquake. We used the Mercator projection with a spatial grid step 0.9 km at latitude 45°N. Computations based on data from DART station 21401 with tsunami duration 30 min were made.

The hindcast was made near the above-mentioned points S-K, K, B and Yu-K. The points are located at distances from the coast, respectively, 1.1, 0.8, 2.1 and 5.4 km. Water depths at these points are 16, 3.2, 30 and 29 m. Locations of these points are shown in Fig. 7.

The results of a hindcast of tsunami waveforms with a duration of 6 - 7 hours made for the above-mentioned points at the Kuril Islands are represented in Fig. 18.

Figure 18.

Tsunami hindcast at the points of the Kuril Islands. The vertical red line – the moment of forecast producing.

The simulated waveform at point Yu-K and the observed tsunami waveform at Yuzhno-Kurilsk are in good agreement. The structure of the tsunami, arrival times of all four waves and wave amplitudes coincided well. According to computations, maximal wave amplitudes should not exceed 1 m; that is in agreement with the observations.

The simulated tsunami waveform at Kurilsk (point K) demonstrates a brightly expressed batch wave structure. The computed wave periods coincide well with recorded periods. According to the forecast, the wave amplitudes should not exceed 0.5 metres.

According to the forecast for Burevestnik (point B), the wave amplitudes should not exceed 1 m. This result agrees with observational data: the wave amplitudes on visual observations were about 1 metre.

At Severo-Kurilsk (point S-K), according to the forecast, increase of amplitude of consequent waves almost up to 1 m is possible. The range of oscillations of a level should not exceed 1.8 - 2 metres. This result agrees well with observational data at Severo-Kurilsk. Based on information from vessels, the swing of sea level was equal to 1.6 m.

Unfortunately, the registration of a tsunami at Burevestnik and Severo-Kurilsk was not conducted. Nevertheless, the results of the computation are confirmed by visual observations and information from vessels in these points.

The time of auxiliary computations for transfer function creation did not exceed 25 min, while the tsunami travel time up to DART station 21401 is equal to 58 min. The auxiliary computations are completed before registration of a tsunami by station 21401. Taking into account the time of tsunami recognition at point 21401 (the first period with a duration of 19 min) the preliminary forecast could be given 78 - 80 min after the main shock of the earthquake. That is 2 hours prior to arrival of the first wave of a tsunami at point S-K, 46-48 minutes prior to arrival of a wave at Yu-K, 40 minutes before tsunami arrival at K and 20 minutes before arrival of the first wave of the tsunami at B. In the latter case the offered hydrophysical method appears to be ineffective.

5. Conclusion

The reciprocity principle, known in acoustics as an asymptotic assumption, is applied to describe non-stationary long ocean waves. The similarity conditions for reciprocal sources, not previously studied, are obtained. It is shown in numerical experiments that, with a quality sufficient for practical applications, the reciprocity principle works quite well, even when asymptotic conditions are not totally satisfied (i.e., when the characteristic time of the wave propagation between the sources is only three to five times greater than the typical wave period or when the characteristic wave length is only two to three times greater than the size of the source area).

The proposed method of computing the tsunami waveform from an unknown source using ocean level data at some points demonstrated its efficiency.

There were questions how the size of auxiliary sources affects the results of computations. Numerical experiments with auxiliary sources of different diameters (with a 20 – 100 km range) gave almost identical results.

There were also questions as to whether it is possible to effectively apply auxiliary axially symmetric sources to create a transfer function. The form of the source influences the directivity of the tsunami. The relief of ocean bottom has a major influence on the evolution of the waves. The use the observational data from the 2009 Simushir, 2010 Chile, 2011 weak Honshu and 2011 Great Honshu tsunamis gives very encouraging results. Despite the fact that the auxiliary source was axially symmetric, the tsunami forecast gives almost identical results based either on the data from reference DART stations located in the same direction as the target sites or from the stations located in other directions.

Despite its approximate character, the proposed method can provide a forecast of tsunami wave parameters for any ocean site with sufficiently high quality and provide real-time operative information for tsunami warning services. The results of this paper demonstrate that the proposed method satisfies the requirements of the warning services and can be effectively used for tsunami warnings.

The numerical experiments described in the paper were conducted for open-ocean target sites. For these sites the method works pretty well, however, the main purpose of TWSs is to provide reliable tsunami forecasts for coastal areas. For this purpose we made auxiliary computations (to evaluate transfer functions) using the detailed bathymetry grid (detailed enough to resolve small-scale resonant features of coastal topography). This allows us to provide the adequate tsunami forecast for the entire coastal zone exposed to arriving tsunami waves. This ensures a straightforward forecast, not only of the leading tsunami waves, but also the following tsunami waves which sometimes are much more intensive than the head waves. The practical interest is to predict tsunami waveforms of long duration, which would take into account secondary waves reflected from island chains, submarine ridges and mainland coasts. Examples of tsunami hindcasts of 2009, 2010, 2011 have shown that the forecast of long duration based on short information from DART stations is possible.

The proposed method of short-term tsunami forecasting uses only seismological information about the location of the earthquake epicentre and works independently of the mechanism of tsunami generation. Numerical experiments show that, despite the anomalous character of the 2011 Great Tohoku tsunami, the hindcast using only level data from DART 21401 without the additional seismological data gives an adequate result. The same example shows that the forecast, based on preliminary computations, is not always effective.

The advantage of the proposed method is that it does not require a pre-computed database of synthetic tsunami waveforms. This is especially important for areas where tsunami early warning systems are just being established and there are no pre-computed databases.

The time of the forecast producing is determined by the time of tsunami propagation up to a point of registration by the DART station and the time of tsunami identification (the first half-cycle, about 15 - 20 minutes). Such a forecast is effective for points for which the difference between the time of tsunami attack and the time of forecast producing is enough for an evacuation of the population and vessels from a threatened zone. For the Kuril Islands this time is equal, conditionally, to 30 minutes.

The comparison of results presented in this paper with some other results based on assimilation of DART data (Japan (East Coast of Honshu) Tsunami, 2011; Tang et al., 2008; Titov, 2009; Wei et al., 2008) shows that the forecast quality obtained by different methods is comparable.

Currently, post-event computations are performed using four non-related software, thence routine work takes a long time. This, unfortunately, does not allow forecasting in real-time.

Appendix A

To verify the condition (9) the problem of wave propagation in a boundless basin is solved. One half of the basin (area A, x < 0) has a constant depth D_A, the other (area B, x > 0) has a constant depth D_B.

Let's the axial–symmetrical initial elevation of a free surface ζ_A0(x-X_A,y-Y_A) with radius R_A is given in area A with centre at a point A. The waveform of a free surface in a point with co-ordinates X_B, Y_B in a right half-space B is sought.

The evolution of waves is described by wave equation with variable water depth (2) with the initial condition ζ(t,x,y)/_t=0 = ζ_A0(x-X_A, y-Y_A) in area S_A and also with conditions on infinity: ζ(t,x,y)→ 0 when |x|, |y|→ ∞.

The Laplace transform on time and Fourier transforms on both co-ordinates x, y are applied to the equation (2) and conditions (3) - (5)

$$f(s,k_x,k_y)={\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}f(t,x,y)\cdot \exp (i{k}_{x}x)\cdot \exp (i{k}_{y}y)\cdot \exp (-st)\cdot dydxdt}}}$$E27

.

So the equation (2) became

$$(\frac{s^2}{c_A^2}+k_x^2+k_y^2)\cdot \zeta (s,k_x,k_y)=\frac{s}{c_A^2}\cdot {\zeta }_{A0}(k_x,k_y),$$E28

where s, k_x, k_y – are parameters of Laplace and Fourier transforms,

c_A – long wave speed on water depth D_A: c_A=(gD_A)^1/2,

$${\zeta }_{A0}(k_x,k_y)=\exp (i{k}_x{X}_A+i{k}_y{Y}_A)\cdot {\displaystyle \underset{S_A}{\int }{\zeta }_{A0}(x^{\prime },y^{\prime })\cdot \exp (i{k}_{x}x\text{'}+i{k}_{y}y\text{'})\cdot d{x}^{\prime }d{y}^{\prime }},$$E29

x', y' – are Cartesian co-ordinates with origin in point A.

The inverse Fourier transform on k_x of function ζ(s,k_x,k_y) from (14):

$${\zeta }_1(s,x,k_y)=\frac{1}{2\pi }\cdot \frac{s}{c_A^2}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{e}^{i{k}_x{X}_A+i{k}_y{Y}_A}{\displaystyle \underset{S_A}{\int }{\zeta }_{A0}(x^{\prime },y^{\prime })\cdot e^{i{k}_x{x}^{\prime }+i{k}_y{y}^{\prime }}d{x}^{\prime }d{y}^{\prime }\frac{e^{-i{k}_{x}x}}{k_x^2+{\chi }_A^2}d{k}_x}},$$E30

where$${\chi }_A^2=\frac{s^2}{c_A^2}+k_y^2$$, is obtained by taking residue. The residue is taken in one pole k_x=-iχ_A, considering the boundedness of the required function.

The result is

$${\zeta }_1(s,x,k_y)=\frac{s}{2c_A^2{\chi }_A}{e}^{i{k}_y{Y}_A}\cdot e^{-{\chi }_A(x-X_A)}\cdot {\displaystyle \underset{}{\int }{\zeta }_{A0}(x^{\prime },y^{\prime })\cdot e^{{\chi }_A{x}^{\prime }+i{k}_y{y}^{\prime }}d{x}^{\prime }d{y}^{\prime }}$$E31

The obtained result is a particular solution of the inhomogeneous equation. By virtue of the fact that the solution is sought in an area with a saltus of depth, it is necessary to take into account a general solution for the homogeneous equation. It is ζ₂(s,x,k_y)=C₂exp(χ_Ax). Only one fundamental solution is taken satisfying the condition ζ₂(s,x,k_y)→0 when x→ -∞.

In a right half-space B a general solution is ζ₃(s,x,k_y)=C₃exp(-χ_Bx) satisfying the condition ζ₃→0 when x→ ∞, where$${\chi }_B^2=\frac{s^2}{c_B^2}+k_y^2$$, c_B=(gD_B)^1/2.

Indeterminate coefficients C₂ and C₃ are sought from the following conditions when x=0

$${\zeta }_1(s,0,k_y)+{\zeta }_2(s,0,k_y)={\zeta }_3(s,0,k_y)$$- continuity of sea level,

$$D_A\cdot \frac{d}{dx}\left({\zeta }_1(s,x,k_y)+{\zeta }_2(s,x,k_y)\right){/}_{x=0}=D_B\cdot \frac{d}{dx}{\zeta }_3(s,x,k_y){/}_{x=0}$$- continuity of stream.

From the two last equations

$$C_3=\frac{2{\chi }_A{D}_A}{{\chi }_A{D}_A+{\chi }_B{D}_B}\cdot \frac{s}{2c_A^2{\chi }_A}\cdot e^{{\chi }_A{X}_A+i{k}_y{Y}_A}\cdot {\displaystyle \underset{S_A}{\int }{\zeta }_{A0}(x^{\prime },y^{\prime })\cdot e^{{\chi }_A{x}^{\prime }+i{k}_y{y}^{\prime }}d{x}^{\prime }d{y}^{\prime }}.$$E32

So, the solution (water level) in the point B(X_B,Y_B) from source A is

$$$$

$${\zeta }_A(s,X_B,k_y)=\frac{s}{g\cdot ({\chi }_A{D}_A+{\chi }_B{D}_B)}\cdot e^{{\chi }_A{X}_A+i{k}_y{Y}_A}\cdot {\displaystyle \underset{S_A}{\int }{\zeta }_{A0}(x^{\prime },y^{\prime })\cdot e^{{\chi }_A{x}^{\prime }+i{k}_y{y}^{\prime }}d{x}^{\prime }d{y}^{\prime }}\cdot e^{-{\chi }_B{X}_B}.$$E33

Inverse Fourier transform on k_y

$${\zeta }_A(s,X_B,Y_B)=\frac{1}{2\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\frac{s}{g\cdot ({\chi }_A{D}_A+{\chi }_B{D}_B)}{e}^{{\chi }_A{X}_A-{\chi }_B{X}_B+i{k}_y{Y}_A}}{\displaystyle \underset{S_A}{\int }{\zeta }_{A0}(x^{\prime },y^{\prime })\cdot e^{{\chi }_A{x}^{\prime }+i{k}_y{y}^{\prime }}d{x}^{\prime }d{y}^{\prime }}\text{}{e}^{-i{k}_y{Y}_B}d{k}_y.$$E34

The regular wave reflection (refraction) from the line of depth changing is hereinafter considered.

It is supposed that |X_A|, |X_B|, |Y_A|, |Y_B|>>R_A. The last integral may be estimated by the steepest descent method (Lavrentiev & Shabat, 1987). The exponent of power Ψ(k_y) = χ_AX_A - χ_BX_B + ik_yY_A - ik_yY_B should be expanded in a power series of k_y near the saddle point. This point is sought from the following condition

$${{\Psi }^{\prime }}_{k_y}(k_y)=\frac{k_y{X}_A}{{\chi }_A}-\frac{k_y{X}_B}{{\chi }_B}+i\cdot (Y_A-Y_B)=0$$E35

.

Taking into account X_A, Y_A<0 appropriate root of this equation is$$k_y=k_{+}=-i\cdot \frac{s}{c_A}\sin {\theta }_A=-i\cdot \frac{s}{c_B}\sin {\theta }_B$$, where θ_A, θ_B – angles of incidence and refraction of a ray, connecting points A and B, in some point on the line x=0.

So, the function Ψ(k_y) and the second derivative Ψ"(k_y) in the saddle point k_y=k₊ are

$$\Psi (k_{+})=\frac{s}{c_A}\cos {\theta }_A\cdot X_A-\frac{s}{c_B}\cos {\theta }_B\cdot X_B+\frac{s}{c_A}\sin {\theta }_A\cdot Y_A-\frac{s}{c_B}\sin {\theta }_B\cdot Y_B=-s\cdot T_{AB},$$E36

where T_AB – propagation time from A to B along a ray and

$${{{\Psi }^{\prime }}^{\prime }}_{k_y}(k_{+})=\frac{c_A{X}_A}{s{\cos }^3{\theta }_A}-\frac{c_B{X}_B}{s{\cos }^3{\theta }_B}.$$E37

Let us assume that ζ_A0(x',y') is an axial symmetric function in circle area S_A with radius R_A. It is convenient to replace the Cartesian co-ordinates on polar (x'=r cosφ, y'=r sinφ) and to enter dimensionless co-ordinates and function (r=r'R_A, ζ_A0(x',y')=a_A ζ'_A0(r') ).

The integrand is substituted by its value in a saddle point and function Ψ(k_y) - by an expansion in a power series near to this point. Thus, the sought-for solution (Laplace transform) is described by expression

$$\begin{array}{l}{\zeta }_A(s,X_B,Y_B)=\\ =\frac{a_A{R}_A^2}{2\pi }\frac{e^{-s\cdot T_{AB}}}{{с}_A\cos {\theta }_A+c_B\cos {\theta }_B}{\displaystyle \underset{0}{\overset{1}{\int }}{{\zeta }^{\prime }}_{A0}(r^{\prime })\cdot {\displaystyle \underset{0}{\overset{2\pi }{\int }}{e}^{\frac{s{R}_A}{c_A}{r}^{\prime }\cos (\varphi -{\theta }_A)}d\varphi }}\cdot r^{\prime }d{r}^{\prime }\cdot \\ \cdot {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\exp (\frac{1}{2}\frac{c_A{X}_A}{s{\cos }^3{\theta }_A}{(k_y+i\frac{s}{c_A}\sin {\theta }_A)}^2-}\frac{1}{2}\frac{c_B{X}_B}{s{\cos }^3{\theta }_B}{(}^{k_y})d{k}_y\end{array}$$E38

If s is imaginary value: s=iσ,- then integral over is $$2\pi \cdot J_0(\frac{\sigma R_A}{c_A}r\text{'})$$ - Bessel function of order 0. So, expression (15) may be rewritten

$$\begin{array}{l}{\zeta }_A(\sigma ,X_B,Y_B)=\\ =a_A{R}_A^2\frac{\exp (-i\sigma \cdot T_{AB})}{{с}_A\cos {\theta }_A+c_B\cos {\theta }_B}{\displaystyle \underset{0}{\overset{1}{\int }}{{\zeta }^{\prime }}_{A0}(r^{\prime })J_0\left(\frac{\sigma R_A}{c_A}r\text{'}\right)}\cdot r^{\prime }d{r}^{\prime }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\exp (-i\frac{1}{2}\frac{c_A{X}_A}{\sigma {\cos }^3{\theta }_A}{(k_y-\frac{\sigma }{c_A}\sin {\theta }_A)}^2}+\\ +i\frac{1}{2}\frac{c_B{X}_B}{\sigma {\cos }^3{\theta }_B}{(}^{k_y})d{k}_y\end{array}$$E39

Analogously, the solution in point A from circle source ζ_B0(x-X_B,y-Y_B) with radius R_B is obtained

$$\begin{array}{l}{\zeta }_B(\sigma ,X_A,Y_A)=\\ =\frac{a_B{R}_B^2}{2\pi }\frac{\exp (-i\sigma \cdot T_{BA})}{{с}_B\cos {\theta }_B+c_A\cos {\theta }_A}{\displaystyle \underset{0}{\overset{1}{\int }}{{\zeta }^{\prime }}_{B0}(r^{\prime })J_0\left(\frac{\sigma R_B}{c_B}r\text{'}\right)}\cdot r^{\prime }d{r}^{\prime }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\exp (i\frac{1}{2}\frac{c_B{X}_B}{\sigma {\cos }^3{\theta }_B}{(k_y+\frac{\sigma }{c_B}\sin {\theta }_B)}^2}-\\ -i\frac{1}{2}\frac{c_A{X}_A}{\sigma {\cos }^3{\theta }_A}{(}^{k_y})d{k}_y\end{array}$$E40

Expressions (16) and (17) are enough for a comparison.

The evaluations of integrals over k_y, obtained with allowance for deformation of an integration contour, are identical.

If the disturbances (sources) are similar, dimensionless function ζ'_A0(r') = ζ'_B0(r').

Then, if volumes of initial elevations are equalled (a_AR²_A = a_BR²_B), functions ζ_A(σ,X_B,Y_B) and ζ_B(σ,X_A,Y_A) will be identical, if the arguments of Bessel functions will be equal:$$\frac{\sigma R_A}{c_A}=\frac{\sigma R_B}{c_B}$$. The similarity condition $$\frac{R_A}{\sqrt{g{D}_A}}=\frac{R_B}{\sqrt{g{D}_B}}$$ (9) follows from here.

Appendix B. Methods of inverse Laplace transform

For a final tsunami evaluation in a specific point A it is necessary to execute the inverse Laplace transforms of (12) or (13). The problems of inverse transform for such expressions are related to ill-posed problems and some regularization methods are required.

The initial functions in (12, 13) are discrete ones, therefore, their Laplace images can be represented by discrete Laplace transforms as a power series with a time step τ, for example, (Doetsch, 1967)

$$\zeta (s,M)=\frac{1-e^{-s\tau }}{s}{\displaystyle \sum _{k=0}^N{\zeta }_k(M)\cdot e^{-ks\tau }}$$E41

,

where ζ_k(M) – sea level in the point M in moments t=kτ, k=0…N, τ – time sample. Nτ is duration of ζ(t,M).

After representing all functions in the form of time series and replacing z=e^-sτ, (13) will be represented as z-transformation

$$\zeta (z,A)={\displaystyle \sum _{k=0}^N{\zeta }_k(A)\cdot z^k}={\displaystyle \sum _{k=0}^{N_1}{\zeta }_k(M)\cdot z^k}\frac{{\displaystyle \sum _{k=0}^{N_2}{\eta }_k(A)\cdot z^k}}{{\displaystyle \sum _{k=0}^{N_3}{\eta }_k(M)\cdot z^k}}.$$E42

The factor$$\frac{1-e^{-s\tau }}{s}$$, being an image of square impulse and being out of interest, hereinafter is omitted.

B.1. The empirical method of inverse Laplace transform: The root-displacement method (Method 1)

After division and multiplication of polynomials in the right part of (18) the polynomial will be derived$$\zeta (z,A)={\displaystyle \sum _{k=0}^N{f}_k{z}^k}={\displaystyle \sum _{k=0}^N{f}_k{e}^{-ks\tau }}$$. The result of inverse Laplace transform of this polynomial is the sequence consisting of δ-functions with factors equal to coefficients of a polynomial. These factors f_k are values of a tsunami level being sought at moments t = kτ, k = 0 … N, after an expected tsunami arrival at point A.

In expression (18) the upper limits of the sums should be identical. Otherwise, in the case of various limits on the result of inverse transforms, there would be doubtful summands.

Unfortunately, the power series obtained in such a way diverges. Two methods are offered below to avoid this effect.

The power series obtained after division and multiplication of polynomials in the right part of (18) is diverging if roots of a polynomial in a denominator have absolute values which are less than 1. The following procedure is applied to avoid a divergence. In a polynomial of a denominator the roots having absolute values less than 1 are substituted by roots with absolute values equalled to 1, but with the same arguments

$${\displaystyle \sum _{n=0}^N{\eta }_n(M)\cdot z^n}\to {\displaystyle \sum _{n=0}^N{\eta }_n(M)\cdot z^n\cdot }\frac{{\displaystyle \prod _{k=0}^K(z-\frac{z_k}{\left|z_k\right|})}}{{\displaystyle \prod _{k=0}^K(z-z_k)}}$$E43

where z_k – root with$$\left|z_k\right|<1$$, K – the number of such roots.

As a result of this procedure, and multiplication and division of power series, the converging power series $${\displaystyle \sum _{n=0}^N{f}_n{z}^n}$$ is obtained, though distorted by a procedure of root-displacement.

For reduction of these distortions, similar operation applied to a denominator is used to obtained series

$${\displaystyle \sum _{n=0}^N{f}_n\cdot z^n}\to {\displaystyle \sum _{n=0}^N{f}_n\cdot z^n}\cdot \frac{{\displaystyle \prod _{k=0}^K\left(z-\frac{z_k}{\left|z_k\right|}\right)}}{{\displaystyle \prod _{k=0}^K\left(z-z_k\right)}}={\displaystyle \sum _{n=0}^N{d}_n\cdot z^n}$$E44

The resulting series $${\displaystyle \sum _{n=0}^N{d}_n{z}^n}$$ is converging and the coefficients d_n are approximate values of the required function.

So the required tsunami waveform in moments t=nτ, n=0…N, is ζ_n(A)=d_n.

B.2. Method of inverse Laplace transform based on the solution of an ill-posed system of linear equations (Method 2)

Representing required function ζ(z,A) as the sum$$\zeta (z,A)={\displaystyle \sum _{k=0}^N{\zeta }_k(A)\cdot z^k}={\displaystyle \sum _{}^{}{\zeta }_k\cdot z^k}$$, the expression (18) can be written as

$${\displaystyle \sum _{k=0}^N{\zeta }_k\cdot z^k}\cdot {\displaystyle \sum _{k=0}^N{\eta }_k(M)\cdot z^k}={\displaystyle \sum _{k=0}^N{\zeta }_k(M)\cdot z^k}\cdot {\displaystyle \sum _{k=0}^N{\eta }_k(A)\cdot z^k},$$E45

where the coefficients ζ_k=ζ_k(A) are required values.

After a multiplication of the sums and the equating of coefficients on identical powers z (in products only N+1 summands are kept) the set of equations can be composed as

$$U\zeta =b,$$E46

where U (U_ij=η_i-j(M)$$i\ge j$$, i=0…N) – trigonal matrix of order N, b – column matrix (b_i=$${\displaystyle \sum _{k=0}^i{\zeta }_k(M)\cdot {\eta }_{i-k}(A)}$$, i=0…N),

ζ - column matrix of unknowns.

Determinant of matrix U in numerical simulations detU=η₀^N+1~10^-200 – 10^-100, when η₀(M) ≈ 0.01 -0.1, N=100. It is obvious that system (19) is ill-posed one (Tikhonov & Arsenin, 1977).

System (19) can be solved using the following method described in (Tikhonov & Arsenin, 1977).

The solution of system (19) is sought from the condition of minimum smoothing functional

$${\Vert U\cdot \zeta -b\Vert }^2+\alpha \cdot {\Vert \zeta \text{}\Vert }^2\to min,\alpha >0$$E47

Norms of matrices are$$\Vert U\Vert =\sqrt{{\displaystyle \sum _{i,j=0}^N{U}_{ij}^2}}$$, $$\Vert \zeta \Vert =\sqrt{{\displaystyle \sum _{i=0}^N{\zeta }_i^2}}$$and$$\Vert b\Vert =\sqrt{{\displaystyle \sum _{i=0}^N{b}_i^2}}$$.

An approximate solution based on the above condition of minimum is a solution of the following system of linear algebraic equations

$$U\text{'}\cdot U\cdot \zeta +\alpha \cdot E\cdot \zeta =U\text{'}\cdot b,$$E48

where U' – transposed matrix, E – unit matrix, – regularization parameter. This regularization parameter may be estimated by the method (Tikhonov & Arsenin, 1977).

B.3. Comparison of the results of two methods of inverse Laplace transform

Figure 19 shows one of the results of application of both methods of the inverse Laplace transform to the computation the tsunami waveform for event of March 9, 2011

Figure 19.

Comparison two methods of inverse Laplace transform.

As seen in Fig. 19, the coincidence of the computed waveforms obtained by different methods (black and red lines), except for the last 5 minutes, is excellent. The real tsunami recording is shown in blue. Further comment on Fig. 19 is not required.

All hindcastings in the present paper were obtained based on Method 2.