# On the contribution of the horizontal sea-bed displacements into the tsunami generation process

Denys Dutykh, Dimitrios Mitsotakis, Leonid Chubarov, Yuriy Shokin

► To cite this version:

Denys Dutykh, Dimitrios Mitsotakis, Leonid Chubarov, Yuriy Shokin. On the contribution of the horizontal sea-bed displacements into the tsunami generation process. Ocean Modelling, 2012, 56, pp.43-56. 10.1016/j.ocemod.2012.07.002. hal-00530999v5

HAL Id: hal-00530999

https://hal.science/hal-00530999v5

Submitted on 11 Jul 2012

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# ON THE CONTRIBUTION OF THE HORIZONTAL SEA-BED DISPLACEMENTS INTO THE TSUNAMI GENERATION PROCESS

DENYS DUTYKH*, DIMITRIOS MITSOTAKIS, LEONID B. CHUBAROV, AND YURIY I. SHOKIN

ABSTRACT. The main reason for the generation of tsunamis is the deformation of the bottom of the ocean caused by an underwater earthquake. Usually, only the vertical bottom motion is taken into account while the horizontal co-seismic displacements are neglected in the absence of landslides. In the present study we propose a methodology based on the well-known Okada solution to reconstruct in more details all components of the bottom coseismic displacements. Then, the sea-bed motion is coupled with a three-dimensional weakly nonlinear water wave solver which allows us to simulate a tsunami wave generation. We pay special attention to the evolution of kinetic and potential energies of the resulting wave while the contribution of the horizontal displacements into wave energy balance is also quantified. Such contribution of horizontal displacements to the tsunami generation has not been discussed before, and it is different from the existing approaches. The methods proposed in this study are illustrated on the July 17, 2006 Java tsunami and some more recent events.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. Mathematical models</td><td>5</td></tr><tr><td>2.1. Dynamic bottom displacements reconstruction</td><td>5</td></tr><tr><td>2.2. The water wave problem with moving bottom</td><td>10</td></tr><tr><td>3. Numerical results</td><td>13</td></tr><tr><td>4. Discussion and conclusions</td><td>22</td></tr><tr><td>Appendix A. Applications to some recent tsunami events</td><td>23</td></tr><tr><td>A.1. Mentawai 2010 tsunami</td><td>23</td></tr><tr><td>A.2. March 11, 2011 Japan earthquake and tsunami</td><td>24</td></tr><tr><td>Acknowledgement</td><td>25</td></tr><tr><td>References</td><td>25</td></tr></table>

## 1. INTRODUCTION

During the years 2004 to 2006 several interesting tsunami events took place in the Indian Ocean. In December 26, 2004, the Great Sumatra-Andaman earthquake ($M_w = 9.1$, cf.

Key words and phrases. tsunami waves; water waves; coseismic displacements; wave energy; finite fault inversion.

<footnote>* Corresponding author. Tel: +33 4 79 75 86 52, Fax: +33 4 79 75 81 42. E-mail: Denys.Dutykh@univ-savoie.fr.</footnote>
C. Ammon et al. (2005), [1] generated the devastating Indian Ocean tsunami refered to in the literature as the Tsunami Boxing Day 2004 (cf. C. Synolakis & E. Bernard (2006), [83]). The local tsunami run-ups from this event exceeded the height of 34 m at Lhoknga in the western Aceh Province. This and other observations led many researchers to ask whether this tsunami was unusualy large for this specific earthquake size. Later it was shown by E.L. Geist et al. (2006) [32] that the Great Sumatra-Andaman earthquake is very similar in terms of local tsunami magnitude to past events of the same size. For example, the 1964 Great Alaska earthquake ($M_w$ = 9.2, cf. H. Kanamori (1970) [49]) demonstrates a similar scaling.

On March 28, 2005 another earthquake occured approximately 110 km to the SE from the Great Sumatra-Andaman earthquake's epicenter. The magnitude of this earthquake was estimated to be $M_w$ = 8.6 ~ 8.7 (cf. [57, 8, 93]). This event triggered a massive evacuation in the surrounding Indian Ocean countries. However, the March 28, 2005 Northern Sumatra earthquake failed to generate a significant tsunami event. The survey teams reported maximum tsunami run-up of 4 m. [32]. This event can be compared with the 1957 Aleutian earthquake ($M_w$ = 8.6, cf. J.M. Johnson and K. Satake (1993) [46]) which produced a maximum tsunami run-up of 15 m (see J.F. Lander (1996) [55]). The deficiency of the March 2005 tsunami is related to the slip concentration in the down-dip part of the rupture zone and to the fact that a substantial part of the vertical displacement occurred in shallow waters or on the substance of the ground [32].

On the other hand, the smaller July 16, 2006 Java earthquake ($M_w$ = 7.8, cf. C.J. Ammon et al. (2006), [2]) generated an unexpectedly destructive tsunami wave which affected over 300 km of south Java coastline and killed more than 600 people, [26]. Field measurements of run-up distributions range uniformly from 5 to 7 m in most inundated areas. However, unexpectedly high run-up values were reported at Nusakambangan Island exceeding the value of 20 m. This tsunami focusing effect could be seemingly ascribed to local site effects and/or to a local submarine landslide/slump. July 16, 2006 Java tsunami can be compared to a similar event occured on June 2, 1994 at the East Java ($M_w$ = 7.6) (see Y. Tsuji et al. (1995), [88]). This 1994 Java tsunami produced more than 200 casualties with local run-up at Rajekwesi slightly exceeding 13 m.

All these examples of recent and historical tsunami events show that there is a big variety of the local tsunami heights and run-up values with respect to the earthquake magnitude $M_w$. It is obvious that the seismic moment $M_0$ of underwater shallow earthquakes is adequate to estimate far-field tsunami amplitudes (see also [68]). However, tsunami wave energy reflects better the local tsunami severity, while the specific run-up distribution depends on bathymetric propagation paths and other site-specific effects. One of the central challenges in the tsunami science is to rapidly assess a local tsunami severity from very first rough earthquake estimations. In the current state of our knowledge false alarms for local tsunamis appear to be unavoidable. The tsunami generation modeling attempts to improve our understanding of tsunami behaviour in the vicinity of the genesis region.

Tsunami generation modeling was initiated in the early sixties by the prominent work of K. Kajiura, [47], who proposed the use of the static vertical sea bed displacement
for the initial condition of the free surface elevation. Classically, the celebrated Okada, [66, 67], and sometimes Mansinha & Smylie<fnref n="1" /> [59, 60] or Gusiakov, [33], solutions are used to compute the coseismic sea bed displacements. This approach is still widely used by the tsunami wave modeling community. However, some progress has been recently made in this direction, [65, 19, 18, 21, 74, 76, 25].

There is a consensus on the importance of the horizontal motion in landslide generated tsunamis, [37, 95, 4, 69]. However, for tsunami waves caused by underwater earthquakes, the horizontal displacements are often disregarded by the tsunami wave modeling community. We can quote here a few publications devoted to tsunami waves such as one by D. Stevenson (2005), [81]:

> "Although horizontal displacements are often larger, they are unimportant for tsunami generation except to the extent that the sloping ocean floor also forces a vertical displacement of the water column."

Or another one by G.A. Ichinose et al. (2000) [41]:

> "The initial lake level values were specified by assuming that the lake surface instantaneously conformed to the vertical displacement of the lake bottom while horizontal velocities were set to zero. The effect of horizontal deformation on the initial condition is neglected here and left for future work."

The authors of this article also tended to neglect the horizontal displacements field in previous tsunami generation studies, [18, 64]. Perhaps, this situation can be ascribed to the work of E. Berg (1970), [7], who showed in the case of the 1964 Alaska's earthquake ($M_w = 9.2$) that the input into the potential energy from the horizontal motion is less than 1.5% of that from the vertical movement.

The attitude to horizontal displacements changed after the prominent work by Y. Tanioka and K. Satake (1996), [84]. According to their computations for the 1994 Java earthquake, the inclusion of horizontal displacements contribution results in 43% increase in maximum initial vertical ground displacement and an increase in the wave amplitudes at the shoreline by 30%. These striking results incited researchers to reconsider the role of the horizontal sea bed motion. Hereafter, various tsunami generation techniques which involve only the vertical displacement field are referred to as the incomplete scenario. On the contrary, the complete tsunami generation takes also into account the horizontal displacements field. The term *complete* in our study refers to a purely kinematic interaction between the moving bottom and the ocean layer in contrast to the horizontal impulse-momentum transfer mechanism [91] applied for the first time to tsunami wave generation in [78].

The approach proposed by Tanioka & Satake (1996) is based on a simple physical consideration: horizontal displacements of a sloping bottom will produce some amount of the vertical motion depending on the slope magnitude. Assuming that the slope is small, this idea can be expressed mathematically by applying once the first order Taylor expansion:

$$u_h := u_1h_x + u_2h_y,$$

<footnote n="1">¹In fact, Mansinha & Smylie solution is a particular case of the more general Okada solution.</footnote>
where $u_{1,2}$ are the horizontal components of the displacements field and $h(x, y)$ is the bathymetry function. The indices $x$ and $y$ denote the spatial partial derivatives. Tanioka & Satake (1996) referred to the quantity $u_h$ as the *vertical displacement of water due to the horizontal movement of the slope*. Then, the quantity $u_h$ is simply added as a correction to the vertical displacements component $u_3(x, y)$. To compute the co-seismic displacements $u_{1,2,3}$ they used the celebrated Okada solution [80, 66, 67] applied to a single rectangular fault segment. In the present study we apply the same solution to multiple sub-faults to achieve higher resolution (see also [24]). More recently, an original approach based on the impulse-momentum principle, [91], has been proposed in [78]. Their genesis theory is based on the idea of the momentum exchange between the slipped continental slope and the water column due to the earthquake. The horizontal bottom motion contributes directly to the ocean kinetic energy as any moving object in water transfers momentum to the fluid. This idea is similar to the classical added mass concept. According to the authors of [78], the continental slope is not an exception from this principle. This implies some special treatment of the water parcels in the vicinity of moving boundaries. For more details we refer to the original publication [78]. The horizontal impulse approach leads to predictions which are somehow different from the present and several previous studies. The methodology described in this article is closer to the ideas of Tanioka & Satake (1996) [84]. Namely, we reconstruct in more details the kinematics of bottom displacements in vertical and horizontal directions, including their evolution in time. We underline that most of tsunami generation studies, including [84], assume the bottom deformation process to be instantaneous, that in most of the cases provides a fairly good approximation. Then, the bottom motion generates some perturbations on the free surface of the water layer. This perturbation form long waves, propagating across the oceans under the gravity force and commonly referred to as *tsunami waves* due to their destructive potential fully realized only in coastal regions, [47, 89, 16, 83].

In the current study we shed some light into the energy transfer process from an underwater earthquake to the implied tsunami wave (in the spirit of the study by D. Dutykh & F. Dias (2009), [20]), while taking into account and quantifying the horizontal displacements contribution into tsunami energy balance. We focus only on the generation stage since the propagation phase and run-up techniques are well understood nowadays at least in the context of Nonlinear Shallow Water equations, [42, 14, 53, 25, 23].

The present study is organized as follows. In Section 2 we present mathematical models used in this study. Specifically, in Sections 2.1 and 2.2 a description of the bottom motion and the fluid layer solution is presented respectively. Some rationale on tsunami wave energy computations is discussed in Section 2.2.1. Numerical results are presented in Section 3. Finally, some important conclusions of this study are outlined in Section 4. In the Appendix of this paper we present the results concerning the tsunami generation of two recent events.
## 2. MATHEMATICAL MODELS

Tsunami waves are caused by a huge and rapid motion of the seafloor due to an underwater earthquake over broad areas in comparison to water depth. There are some other mechanisms of tsunami genesis such as underwater landslides, for example. However, in this study we focus on the purely seismic mechanism which occurs most frequently in nature.

Hydrodynamics and seismology are only weakly coupled in the tsunami generation problem. Namely, the released seismic energy is partially transmitted to the ocean through the sea bed deformation while the ocean has obviously no influence onto the rupturing process. Consequently, our problem is reduced to two relatively distinct questions:

(1) Reconstruct the sea bed deformation $h = h(\vec{x}, t)$ caused by the seismic event under consideration

(2) Compute the resulting free surface motion

The answers to these questions are analyzed in Sections 2.1 and 2.2 respectively.

### 2.1. Dynamic bottom displacements reconstruction.
Traditionally, the free surface initial condition for various tsunami propagation codes (see [85, 42, 43]) is assumed to be identical to the static vertical deformation of the ocean bottom, [48]. This assumption is classically justified by the three following reasons:

(1) Tsunamis are long waves. In this regime the vertical acceleration is neglected with respect to the gravity force

(2) The sea bed deformation is assumed to be instantaneous. It is based on the comparison of gravity wave speed (200 m/s for water depth of 4 km) and the seismic wave celerity (≈ 3600 m/s)

(3) The effect of horizontal bottom motion is negligible for tsunami generation since the bathymetry has in general mild slope (≈ 10%), [7]

It is worth to note that nonhydrostatic effects as well as finite time source duration have been modeled in several recent studies [87, 19, 28, 52, 21, 29]. Some attempts have also been made to take into account the horizontal displacements contribution, [84, 32, 78].

In the present study we relax all three classical assumptions. The bottom deformation is reconstructed using the finite fault solution as it was suggested in our previous study [24]. However, we extend the previous construction to take into account the horizontal displacements contribution as well. The finite fault solution is based on the multi-fault representation of the rupture, [5, 45]. The rupture complexity is reconstructed using the inversion of seismic data. Fault's surface is parametrized by multiple segments with variable local slip, rake angle, rise time and rupture velocity. The inversion is performed in an appropriate wavelet transform space. The objective function is a weighted sum of $L_1, L_2$ norms and some correlative functions. With this approach seismologists are able to recover rupture slip details, [5, 45]. This available seismic information is exploited hereafter to compute the sea bed displacements produced by an underwater earthquake with better geophysical resolution. Several other multiple segments sources have been used in [31, 94, 61].
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Surface projection of the fault's plane and the ETOPO1 bathymetric map of the region under consideration. The symbol ⋆ indicates the hypocenter's location at (107.345°, -9.295°). The local Cartesian coordinate system is centered at the point (108°, -10°). This region is located between (106°, -8°) and (110°, -12°). The colorbar indicates the water depth in meters below the still water level (z = 0). In the region under consideration the depth varies from 20 to 7100 meters.</figcaption></figure>

**Remark 1.** In a few studies an attempt has been made to reconstruct the seismic source from tsunami tide gauge records, [71, 72, 97, 30, 40, 75]. This approach seems to be very promising and in future a joint combination of seismic and hydrodynamic inversions should be used for the successful reconstruction of appropriate initial conditions.

Let us describe the way of how the sea bed displacements are reconstructed. In this reconstruction procedure we follow the great lines of our previous study [24], while adding new ingredients concerning the horizontal displacements field reconstruction. We illustrate the proposed approach in the case of the July 17, 2006 Java tsunami for which the finite fault solution is available, cf. [44, 70]. It was a relatively slow earthquake and thus, atypical. However, we assume that the slow rupturing process was well resolved by the finite fault inversion algorithm.

The fault is considered to be the rectangle with vertices located at (109.20508° (Lon), -10.37387° (Lat), 6.24795 km (Depth)), (106.50434°, -9.45925°, 6.24795 km), (106.72382°, -8.82807°, 19.79951 km), (109.42455°, -9.74269°, 19.79951 km) (see Figure 1). The fault's plane is conventionally divided into $N_x = 21$ subfaults along strike and $N_y = 7$ subfaults
<table><caption>TABLE 1. Geophysical parameters used to model elastic properties of the subduction zone in the region of Java.</caption><tr><td>P-wave celerity $c_p$, m/s</td><td>6000</td></tr><tr><td>S-wave celerity $c_s$, m/s</td><td>3400</td></tr><tr><td>Crust density $\rho$, kg/m³</td><td>2700</td></tr><tr><td>Dip angle, $\delta$</td><td>10.35°</td></tr><tr><td>Slip angle (CW from N)</td><td>288.94°</td></tr></table>

down the dip angle, leading to the total number of $N_x \times N_y = 147$ equal segments. Parameters such as the subfault location $(x_c, y_c)$, the depth $d_i$, the slip $u_i$ and the rake angle $\phi_i$ for each segment can be found in [44] (see also Appendix II, [24]). The elastic constants and parameters such as dip and slip angles, which are common to all subfaults, are given in Table 1. We underline that the slip angle is measured conventionally in the counter-clockwise direction from the North. The relations between the elastic wave celerities $c_p, c_s$ and Lamé coefficients $\lambda, \mu$ used in Okada's solution are classical and can also be found in Appendix III, [24].

One of the main ingredients in our construction is the so-called Okada solution, [66, 67], which is used in the case of an active fault of small or intermediate size. The success of this solution may be ascribed to the closed-form analytical expressions which can be effectively used for various modeling purposes involving co-seismic deformation.

**Remark 2.** The celebrated Okada solution, [66, 67], is based on two main ingredients — the dislocation theory of Volterra [92, 80] and Mindlin's fundamental solution for an elastic half-space, [63]. Particular cases of this solution were known before Okada's work, for example the well-known Mansinha and Smylie's solution, [59, 60]. Usually all these particular cases differ by the choice of the dislocation and the Burger's vector orientation, [73]. We recall the basic assumptions behind this solution:

* Fault is immersed into the linear homogeneous and isotropic half-space

* Fault is a Volterra's type dislocation

* Dislocation has a rectangular shape

For more information on Okada's solution we refer to [19, 16, 18] and the references therein.

The trace of the Okada's solution at the sea bottom (substituting $z = 0$ in the geophysical coordinate system) for each subfault will be denoted $O_i^{(j)}(\vec{x}; \delta, \lambda, \mu, ...)$, where $\delta$ is the dip angle, $\lambda, \mu$ are the Lamé coefficients and dots denote the dependence of the function $O_i^{(j)}(\vec{x})$ on other 8 parameters, cf. [19]. The index $i$ takes values from 1 to $N_x \times N_y$ and denotes the corresponding subfault segment, while the superscript $j$ is equal to 1 or 2 for the horizontal displacements and to 3 for the vertical component. Hereafter we will adopt the short-hand notation $O_i^{(j)}(\vec{x})$ for the $j$th displacement component of the Okada's solution for the $i$th segment having in mind its dependence on various parameters.

Taking into account the dynamic characteristics of the rupturing process, we make some further assumptions on the time dependence of the displacement fields. The finite fault
solution provides us with two additional parameters concerning the rupture kinematics — the rupture velocity $v_r$ and the rise time $t_r$ which are equal to 1.1 km/s and 8 s for July 17, 2006 Java event respectively. The epicenter is located at the point $\vec{x}_e = (107.345^\circ, -9.295^\circ)$, cf. [44]. Given the origin $\vec{x}_e$, the rupture velocity $v_r$ and $i^{th}$ subfault location $\vec{x}_i$, we define the subfault activation times $t_i$ needed for the rupture to achieve the corresponding segment $i$ by the formulas:

$$t_i = \frac{\|\vec{x}_e - \vec{x}_i\|}{v_r}, \quad i = 1, \dots, N_x \times N_y.$$

Consequently, for the sake of simplicity in our study we assume the rupture propagation along the fault to be isotropic and homogeneous. However, some more sophisticated approaches including possible heterogeneities exist (see e.g. [13]). We will also follow the pioneering idea of J. Hammack, [36, 35], developed later in [87, 86, 19, 22, 52], where the maximum bottom deformation is achieved during some finite time (the so-called *rise time*) according to an appropriately chosen dynamic scenario. Various scenarios used in practice (instantaneous, linear, trigonometric, exponential, etc) can be found in [36, 50, 35, 22, 19]. In this study we adopt the trigonometric scenario which is given by the following formula:

$$T(t) = \mathcal{H}(t - t_r) + \frac{1}{2}\mathcal{H}(t)\mathcal{H}(t_r - t)(1 - \cos(\pi t / t_r)), \quad (2.1)$$

where $\mathcal{H}(t)$ is the Heaviside step function. This scenario has the advantage to have also the first derivative continuous at the activation time $t = 0$. However, for comparative purposes sometimes we will use also the so-called exponential scenario [50]:

$$T_e(t) = \mathcal{H}(t)(1 - e^{-\alpha t}), \quad \alpha := \frac{\log(3)}{t_r}.$$

For illustrative purposes both dynamic scenarios are represented in Figure 2.

We sum up together all the ingredients proposed above to reconstruct dynamic displacements field $\vec{u} = (u_1, u_2, u_3)$ at the sea bottom:

$$u_j(\vec{x}, t) = \sum_{i=1}^{N_x \times N_y} T(t - t_i)\mathcal{O}_i^{(j)}(\vec{x}).$$

**Remark 3.** We would like to underline here the asymptotic behaviour of the sea bed displacements. By definition of the trigonometric scenario (2.1) we have $\lim_{t \to +\infty} T(t) = 1$. Consequently, the sea bed deformation will attain fast its state which consists of the linear superposition of subfaults contributions:

$$u_j(\vec{x}, t) = \sum_{i=1}^{N_x \times N_y} \mathcal{O}_i^{(j)}(\vec{x}).$$

Finally, we can predict the sea bed motion by taking into account horizontal and vertical displacements:

$$h(\vec{x}, t) = h_0(\vec{x} - \vec{u}_{1,2}(\vec{x}, t)) - u_3(\vec{x}, t), \quad (2.2)$$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Trigonometric and exponential dynamic scenarios for $t_r = 1$ s (see J. Hammack (1973), [35]).</figcaption></figure>

where $h_0(\vec{x})$ is a function which interpolates<fnref n="2" /> the static bathymetry profile given e.g. by the ETOPO1 database (see Figures 1 and 4).

**Remark 4.** In some studies where horizontal displacements were taken into account (cf. [84, 3, 78, 6]), the first order Taylor expansion was permanently applied to the bathymetry representation formula (2.2) to give:

$$h(\vec{x}, t) \approx h_0(\vec{x}) - \vec{u}_{1,2}(\vec{x}, t) \cdot \nabla_{\vec{x}} h_0(\vec{x}) - u_3(\vec{x}, t) \quad (2.3)$$

We prefer not to follow this approximation and to use the exact formula (2.2) since it is valid for all slopes (see Figure 4 for Java bathymetry). Another difficulty lies in the estimation of the bathymetry gradient $\nabla_{\vec{x}}h_0(x)$ required by Taylor's formula (2.3). The application of any finite differences scheme to the measured $h_0(\vec{x})$ leads to an ill-posed problem. Consequently, one needs to apply extensive smoothing procedures to the raw data $h_0(\vec{x})$ which induces an additional loss in accuracy.

In the present study we do not completely avoid the computation of the static bathymetry gradient $\nabla_{\vec{x}}h_0(\vec{x})$ since the kinematic bottom boundary condition (2.7) involves the time derivative of the bathymetry function:

$$\partial_t h = -\nabla_{\vec{x}}h_0(\vec{x}) \cdot \partial_t \vec{u}_{1,2}(\vec{x}, t) - \partial_t u_3(\vec{x}, t).$$

<footnote n="2">²In our numerical simulations presented below we use the MATLAB TriScatteredInterp class to interpolate the static bathymetry values given by ETOPO1 database.</footnote>
However, the last formula is exact and it is obtained by a straightforward application of the chain differentiation rule.

Another possibility could be to consider static horizontal displacements $\vec{u}_{1,2}(\vec{x})$ thus keeping dynamics only in the vertical component $u_3(\vec{x}, t)$. However, we do not choose this option in this work.

In the next section we will present our approach in coupling this dynamic deformation with the hydrodynamic problem to predict waves induced on the free surface of the ocean.

**2.2. The water wave problem with moving bottom.** We consider the incompressible flow of an ideal fluid with constant density $\rho$ in the domain $\Omega \subseteq \mathbb{R}^2$. The horizontal independent variables will be denoted by $\vec{x} = (x, y)$ and the vertical one by $z$. The origin of the Cartesian coordinate system is traditionally chosen such that the surface $z = 0$ corresponds to the still water level. The fluid domain is bounded below by the bottom $z = -h(\vec{x}, t)$ and above by the free surface $z = \eta(\vec{x}, t)$. Usually we assume that the total depth $H(\vec{x}, t) := h(\vec{x}, t) + \eta(\vec{x}, t)$ remains positive $H(\vec{x}, t) \ge h_0 > 0$ under the system dynamics $\forall t \in [0, T]$. The sketch of the physical domain is shown in Figure 3.

**Remark 5.** Classically in water wave modeling, we make the assumption that the free surface is a graph $z = \eta(\vec{x}, t)$ of a single-valued function. It means in practice that we exclude some interesting phenomena, (e.g. wave breaking phenomena) which are out of the scope of this modeling paradigm.

The governing equations of the classical water wave problem are the following, [54, 82, 62, 96]:

$$ \Delta\phi = \nabla^2\phi + \partial_{zz}^2\phi = 0, \quad (\vec{x}, z) \in \Omega \times [-h, h], \quad (2.4) $$

$$ \partial_t\eta + \nabla\phi \cdot \nabla\eta - \partial_z\phi = 0, \quad z = \eta(\vec{x}, t), \quad (2.5) $$

$$ \partial_t\phi + \frac{1}{2}|\nabla\phi|^2 + \frac{1}{2}(\partial_z\phi)^2 + g\eta = 0, \quad z = \eta(\vec{x}, t), \quad (2.6) $$

$$ \partial_t h + \nabla\phi \cdot \nabla h + \partial_z\phi = 0, \quad z = -h(\vec{x}, t), \quad (2.7) $$

where $\phi$ is the velocity potential, $g$ the acceleration due to gravity force and $\nabla = (\partial_x, \partial_y)$ denotes the gradient operator in horizontal Cartesian coordinates. The fluid incompressibility and flow irrotationality assumptions lead to the Laplace equation (2.4) for the velocity potential $\phi(\vec{x}, z, t)$.

The main difficulty of the water wave problem lies on the boundary conditions. Equations (2.5) and (2.7) express the free surface and bottom impermeability, while the Bernoulli condition (2.6) expresses the free surface isobarity respectively.

Function $h(\vec{x}, t)$ represents the ocean's bathymetry (depth below the still water level, see Figure 3) and is assumed to be known. The dependence on time is included in order to take into account the bottom motion during an underwater earthquake [16, 19, 22, 52, 18].

For the exposition below we will need also to compute unitary exterior normals to the fluid domain. The normals at the free surface $S_f$ and at the bottom $S_b$ are given by the
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Sketch of the physical domain.</figcaption></figure>

following expressions respectively:

$$ \hat{n}_f = \frac{1}{\sqrt{1 + |\nabla \eta|^2}} \begin{vmatrix} -\nabla \eta & 1 \end{vmatrix}, \quad \hat{n}_b = \frac{1}{\sqrt{1 + |\nabla h|^2}} \begin{vmatrix} -\nabla h & -1 \end{vmatrix}. $$

In 1968 V. Zakharov proposed a different formulation of the water wave problem based on the trace of the velocity potential at the free surface [98]:

$$ \varphi(\vec{x}, t) := \phi(\vec{x}, \eta(\vec{x}, t), t). $$

This variable plays a role of the generalized momentum in the Hamiltonian description of water waves [98, 15]. The second canonical variable is the free surface elevation $\eta$.

Another important ingredient is the normal velocity at the free surface $v_n$ which is defined as:

$$ v_n(\vec{x}, t) := \sqrt{1 + |\nabla \eta|^2} \frac{\partial \phi}{\partial \hat{n}_f} \bigg|_{z=\eta} = (\partial_z \phi - \nabla \phi \cdot \nabla \eta)|_{z=\eta}. \quad (2.8) $$

Kinematic and dynamic boundary conditions (2.5), (2.6) at the free surface can be rewritten in terms of $\varphi$, $v_n$ and $\eta$ [12, 11, 27]:

$$ \begin{aligned} \partial_t \eta - \mathcal{D}_\eta(\varphi) &= 0, \\ \partial_t \varphi + \frac{1}{2} |\nabla \varphi|^2 + g\eta - \frac{1}{2(1+|\nabla \eta|^2)} \left[ \mathcal{D}_\eta(\varphi) + \nabla \varphi \cdot \nabla \eta \right]^2 &= 0. \end{aligned} \quad (2.9) $$
Here we introduced the so-called Dirichlet-to-Neumann operator (D2N) [10, 11] which maps the velocity potential at the free surface $\varphi$ to the normal velocity $v_n$:

$$ \mathcal{D}_{\eta}: \quad \varphi \mapsto v_n = \sqrt{1 + |\nabla \eta|^2} \frac{\partial \phi}{\partial \hat{n}_f} \Big|_{z=\eta} \\ \left| \begin{array}{l} \nabla^2 \phi + \partial_{zz}^2 \phi = 0, \quad (\vec{x}, z) \in \Omega \times [-h, \eta], \\ \dot{\phi} = \varphi, \quad z = \eta, \\ \sqrt{1 + |\nabla h|^2} \frac{\partial \phi}{\partial \hat{n}_b} = \partial_t h, \quad z = -h. \end{array} \right. $$

The name of this operator comes from the fact that it makes a correspondance between Dirichlet data $\varphi$ and Neumann data $\sqrt{1 + |\nabla \eta|^2} \frac{\partial \phi}{\partial \hat{n}_f} \Big|_{z=\eta}$ at the free surface.

So, the water wave problem can be reduced to a system of two PDEs (2.9) governing the evolution of the canonical variables $\eta$ and $\varphi$. For the tsunami generation problem we approximate and we compute efficiently the D2N map $\mathcal{D}_{\eta}(\varphi)$ using the Weakly Nonlinear (WN) model described in [24]. This relies on the approximate solution of the 3D Laplace equation in a perturbed strip-like domain using the Fourier transform ($\hat{\varphi} := \mathcal{F}[\varphi]$, $\eta = \mathcal{F}^{-1}[\eta]$):

$$ \hat{\mathcal{D}}_{\eta}(\varphi) = \hat{\varphi}|_{\vec{k}} \tanh(|\vec{k}|H) + \hat{f} \operatorname{sech}(|\vec{k}|H) - \mathcal{F}[\mathcal{F}^{-1}[i\vec{k}\hat{\varphi}] \cdot \mathcal{F}^{-1}[i\vec{k}\hat{\eta}]] $$

where $\vec{k}$ is the wavenumber and $f$ is the bathymetry forcing term

$$ f(\vec{x}, t) := -\partial_t h - \nabla \phi|_{z=-h} \cdot \nabla h. $$

Several details on the time integration procedure can be also found in [24]. The resulting method is only weakly nonlinear and analogous at some point to the first order approximation model proposed in [34]. As the hydrodynamical models, Tanioka and Satake (1996) [84] used the Linear Shallow Water equations [77], while Song et al. (2008) employed a hydrostatic ocean circulation model [79].

Recently, it was shown that a tsunami generation process is essentially linear, [52, 76]. However, our WN approach will take into account some nonlinear effects when they become important, for example when rapid changes in the bathymetry are present. This is possible when the generation region involves a wide range of depths from deep to shallow regions (see Figures 1 and 4).

### 2.2.1. *Tsunami wave energy.*
In this study we are more particularly interested in the evolution of the generated wave energy [20]. In the case of free surface incompressible flows, the kinetic and potential energies, denoted by $K$ and $\Pi$ respectively, are completely determined by the velocity field and the free surface elevation:

$$ K(t) := \frac{\rho}{2} \int_{-h}^{\eta} \iint_{\Omega} |\nabla \phi|^2 d\vec{x} dz, \quad \Pi(t) := \frac{\rho g}{2} \iint_{\Omega} \eta^2 d\vec{x}. $$

The definition of the kinetic energy $K(t)$ involves an integral over the three dimensional physical domain $\Omega \times [-h, \eta]$. We can reduce the integral dimension using the fact that the
velocity potential $\phi$ is a harmonic function:

$$|\nabla \phi|^2 = \nabla \cdot (\phi \nabla \phi) - \underbrace{\phi \Delta \phi}_{=0} \equiv \nabla \cdot (\phi \nabla \phi).$$

Consequently, the kinetic energy can be rewritten as follows:

$$K(t) = \frac{\rho}{2} \iint_{S_f+S_b} \phi \nabla \phi \cdot \hat{n} \, d\sigma = \underbrace{\frac{\rho}{2} \iint_{\Omega} \varphi \mathcal{D}_{\eta}(\varphi) \, d\vec{x}}_{(I)} + \underbrace{\frac{\rho}{2} \iint_{\Omega} \check{\phi} \partial_t h \, d\vec{x}}_{(II)}$$

where $\check{\phi}$ denotes the trace of the velocity potential at the bottom $\phi|_{z=-h}$ (see D. Clamond & D. Dutykh (2012), [9]). In order to obtain the last equality we used the free surface and the bottom kinematic boundary conditions (2.6), (2.7). The first integral (I) is classical and represents the change of kinetic energy under the free surface motion while the second one (II) is the forcing term due to the bottom deformation. The total energy<fnref n="3" /> is defined as the sum of kinetic and potential ones:

$$E(t) := K(t) + \Pi(t) = \frac{\rho}{2} \iint_{\Omega} \varphi \mathcal{D}_{\eta}(\varphi) \, d\vec{x} + \frac{\rho}{2} \iint_{\Omega} \check{\phi} \partial_t h \, d\vec{x} + \frac{\rho g}{2} \iint_{\Omega} \eta^2 \, d\vec{x}.$$

Below we will compute the evolution of the kinetic, potential and total energies beneath moving bottom.

### 3. NUMERICAL RESULTS

The proposed approach will be directly illustrated on the Java 2006 event. The July 17, 2006 Java earthquake involved thrust faulting in the Java's trench and generated a tsunami wave that inundated the southern coast of Java, [2, 26]. The estimates of the size of the earthquake, [2], indicate a seismic moment of $6.7 \times 10^{20} \text{ N} \cdot \text{m}$, which corresponds to the magnitude $M_w = 7.8$. Later this estimation was refined to $M_w = 7.7$, [44]. Like other events in this region, Java's event had an unusually low rupture speed of 1.0 – 1.5 km/s (we take the value of 1.1 km/s according to the finite fault solution [44]), and occurred near the up-dip edge of the subduction zone thrust fault. According to C. Ammon et al, [2], most aftershocks involved normal faulting. The rupture propagated approximately 200 km along the trench with an overall duration of approximately 185 s. The fault's surface projection along with ocean ETOPO1 bathymetric map are shown in Figures 1 and 4. We note that Indian Ocean's depth of the region considered in this study varies between 7186 and 20 meters in the shallowest regions which may imply local importance of nonlinear effects.

**Remark 6.** We have to mention that the finite fault inversion for this earthquake was also performed by the Caltech team, [70]. The main differences with the USGS inversion consist on the employed dataset. To our knowledge, A. Ozgun Konca and his collaborators include also displacements measured with GPS-based techniques. Consequently, they came

<footnote n="3">³We note that the total energy is not conserved during the tsunami generation phase due to the forcing term (II) coming from the bottom kinematic boundary condition.</footnote>
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Side view of the bathymetry, cf. also Figure 1.</figcaption></figure>

<table><caption>TABLE 2. Values of physical and numerical parameters used in simulations.</caption><tr><td>Ocean water density, ρ, kg/m<sup>3</sup></td><td>1027.0</td></tr><tr><td>Gravity acceleration, g, m/s<sup>2</sup></td><td>9.81</td></tr><tr><td>Epicenter location (Lon, Lat)</td><td>(107.345°, -9.295°)</td></tr><tr><td>Rupture velocity, v<sub>r</sub>, km/s</td><td>1.1</td></tr><tr><td>Rise time, t<sub>0</sub>, s</td><td>8.0</td></tr><tr><td>Number of Fourier modes in x</td><td>256</td></tr><tr><td>Number of Fourier modes in y</td><td>256</td></tr></table>

to the conclusion that the July 17, Southern Java earthquake magnitude was $M_w = 7.9$. The energy of a tsunami wave generated according to this solution will be discussed below.

The numerical solutions presented below are given by the Weakly Nonlinear (WN) model. A uniform grid of 256 × 256 points<fnref n="4" /> is used in all computations below. The time step $\Delta t$ is chosen adaptively according to the RK(4,5) method proposed in [17]. The problem is integrated numerically during $T = 255$ s which is a sufficient time interval for the bottom to take its final shape (< 220 s) and of the resulting tsunami wave to enter into

<footnote n="4">⁴Since we use a pseudo-spectral method, the convergence is expected to be exponential and this number of harmonics should be sufficient to capture all scales important for phenomena that we consider here.</footnote>
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Contribution of the horizontal displacements field into the sea bed deformation expressed in percentage of the maximum vertical displacement. The grey color corresponds to negligible values while white and black zones show the most important contributions. The colorbar shows the percentage of the signed difference going from -21% to +7%.</figcaption></figure>

the propagation stage. The values of various physical and numerical parameters used in simulations are given in Table 2.

We begin our numerical investigations by quantifying the contribution of horizontal displacements into the sea bed deformation process. For this purpose we consider the difference $d_h(\vec{x})$ between the deformed bottom under the action of only horizontal displacements in their steady state ($t \to +\infty$) and the initial configuration:

$$d_h(\vec{x}) := \frac{h_0(\vec{x} - \vec{u}_{1,2}(\vec{x}, t)|_{t \to +\infty}) - h_0(\vec{x})}{\max_{\vec{x}, t \to +\infty} |u_3(\vec{x}, t)|}.$$

In Figure 5 we present the quantitative effect of horizontal displacements relative to the maximum vertical displacement $\max_{\vec{x}} |u_3(\vec{x}, t = +\infty)|$. The computations we performed show that the maximum amplitude of the bottom variation due to the action of horizontal displacements reaches 21% of the maximum amplitude of the vertical displacement. In practice it means that locally (depending on the bathymetry shape and the slip distribution) we cannot completely neglect the effect of horizontal motion.

The next step consists in quantifying the impact of horizontal displacements onto free surface motion. We put six numerical wave gauges at the following locations: (a) (107.2°, -9.388°), (b) (107.4°, -9.205°), (c) (107.6°, -9.648°), (d) (107.7°, -9.411°), (e) (108.3°, -10.02°), (f) (108.2°, -9.75°). The locations of the wave gauges are represented on Figure
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Location of the six numerical wave gauges (indicated by the symbol ◊) superposed with the steady state coseismic bottom displacement (only the vertical component in meters is represented here).</figcaption></figure>

6 along with the static sea bed displacement. Wave gauges are intentionally put in places where the largest waves are expected. Synthetic wave gauge records are presented in Figure 7. We consider the following four scenarios:

• trigonometric, complete (blue solid line)

• trigonometric, incomplete (black dashed line)

• exponential, complete (blue dash-dotted line)

• exponential, incomplete (black dotted line)

The importance of the nonlinear effects have already been addressed in several previous studies, [52, 76]. In the framework of the Weakly Nonlinear approach we studied this question in a recent companion paper [24] using only the vertical deformation. A fairly good agreement has been observed with the Cauchy-Poisson (linear, fully dispersive) formulation in accordance with preceding results, cf. [52, 76]. Consequently, in the present study we decided to focus on the further comparison between complete/incomplete approaches and exponential/trigonometric scenarios, [35, 22, 64]. These results are discussed hereafter.

We note that the trigonometric scenario leads in general to slightly larger amplitudes than the exponential bottom motion. It is not surprising since the exponential scenario prescribes smoother and less rapid change in the bottom for the same rise time parameter value (see Figure 2). Later we will consider the trigonometric scenario unless otherwise noted.

Then, it can be seen that in most cases horizontal displacements lead to an increase in the wave amplitude but not always as it can be observed in Figure 7(f) (the last wave gauge located at (108.2°, -9.75°)). We investigated more thoroughly this question. Figure 8 shows the relative difference between free surface elevations computed according to
<figure><img src="image_8.png" /><figcaption>(a) Gauge at (107.2°, -9.388°)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Gauge at (107.4°, -9.205°)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(c) Gauge at (107.6°, -9.648°)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(d) Gauge at (107.7°, -9.411°)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(e) Gauge at (108.3°, -10.02°)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(f) Gauge at (108.2°, -9.75°)</figcaption></figure>

FIGURE 7. Free surface elevations computed numerically at six wave gauges located approximately in local extrema of the static bottom displacement. The vertical axis is represented in meters and time is given in seconds. The black lines correspond to the wave generated only by the vertical displacements (incomplete generation) while blue lines take also into account the horizontal displacements contribution (complete scenario).
<figure><img src="image_14.png" /><figcaption>FIGURE 8. Relative difference between the free surface elevation at $t_e = 220$ seconds computed according to the complete $\eta_c(\vec{x}, t_e)$ (with horizontal displacements) and incomplete $\eta_i(\vec{x}, t_e)$ (only vertical component) tsunami generation scenarios. The bottom moves according to the trigonometric scenario (2.1). The vertical scale is given in percents of the maximum amplitude of the incomplete scenario — $\max_{\vec{x}} \eta_i(\vec{x}, t_e)$.</figcaption></figure>

complete $\eta_c(\vec{x}, t_e)$ and incomplete $\eta_i(\vec{x}, t_e)$ scenarios. More precisely, we plot the following quantity (for the trigonometric scenario):

$$d(\vec{x}) := \frac{\eta_c(\vec{x}, t_e) - \eta_i(\vec{x}, t_e)}{\max_{\vec{x}} |\eta_i(\vec{x}, t_e)|}, \quad t_e = 220 \text{ s.}$$

Time $t_e = 220$ s has been chosen because at that moment the bottom has been stabilized and the waves enter into the free propagation regime. In Figure 8 the grey color corresponds to the zero value of the difference $d(\vec{x})$, while the white color shows regions where the wave is amplified by horizontal displacements by approximately 10%. On the contrary, black zones show an attenuation effect of horizontal sea bed motion (about -5%). Recall that all values are given in terms of the maximum amplitude $\max_{\vec{x}} |\eta_i(\vec{x}, t_e)|$ percentage of the incomplete generation approach. Some connection with the results presented in Figure 5 can be noticed.

Finally, we study the evolution of kinetic, potential and total energies during the tsunami generation process described in Section 2.2.1. Specifically we are interested in quantifying the contribution of horizontal displacements into tsunami energy balance. Figure 9 shows the evolution of potential (9(a)) and kinetic (9(b)) energies for four cases already mentioned
<figure><img src="image_15.png" /><figcaption>FIGURE 9. Energy evolution during our simulations in the complete (blue solid line) and incomplete (black dotted line) scenarios. Note that scales are different on the left and right images. The time *t* is given in seconds.</figcaption></figure>

above. Here again blue lines refer to complete generation scenarios while black lines – to vertical displacements only. Consecutive peaks in the kinetic energy come from the activation of new fault segments in accordance with the rupture propagation along the fault. The energy curves vary in an analogous way with the tide gauges records (see Figure 7). Namely, the trigonometric scenario leads to slightly higher energies than the exponential one. However, once the bottom motion stops, this difference becomes negligible. The inclusion of horizontal displacements has a much more visible effect with higher energies.

We performed also the same simulation but using the finite fault solution obtained by the Caltech team, [70], and the trigonometric scenario. The results concerning the kinetic and potential energies are presented in Figure 10. It is interesting to note that the potential energy evolution in the Caltech scenario appears to be smoother especially after $t = 150$ s. However, the peaks corresponding to subfaults activation times are equally present in kinetic energies. The magnitudes of kinetic and potential energies predicted according to the Caltech inversion are approximatively 10 times higher than corresponding USGS results.

We have to underline that inversions performed by the two different finite fault algorithms lead to very different results. Our method is operational with both of them. It is particularly interesting to compare the total energies evolution predicted by USGS and Caltech inversions. These results are presented in Figure 11. The Caltech version of the rupturing process generates a tsunami wave with much higher energy (computed after the end of the bottom motion when a tsunami enters in the propagation regime):

* USGS, complete trigonometric: $E_c = 2.16 \times 10^{12}$ J,

* USGS, complete exponential: $E_c = 2.11 \times 10^{12}$ J,

* USGS, incomplete trigonometric: $E_i = 1.95 \times 10^{12}$ J,
<figure><img src="image_16.png" /><figcaption>FIGURE 10. Energy evolution during our simulations in the complete (blue solid line) and incomplete (black dotted line) scenarios using the finite fault inversion by Caltech Tectonics Observatory [70]. Note that the scales are different on the left and right images. The time *t* is given in seconds.</figcaption></figure>

• USGS, incomplete exponential: $E_i = 1.90 \times 10^{12}$ J,

• Caltech, complete trigonometric: $E_c = 2.28 \times 10^{13}$ J,

• Caltech, incomplete trigonometric: $E_c = 1.83 \times 10^{13}$ J.

This result is not surprising since according to the USGS solution the Java 2006 event magnitude is $M_w = 7.7$ and according to the Caltech team, $M_w = 7.9$, which means a huge difference since the scale is logarithmic.

**Remark 7.** We would like to note an interesting property. For linear waves it can be rigorously shown the exact equipartition property between the kinetic and potential energies [96]. When the nonlinearities are included, this equidistribution property is only approximate. One can observe that at the final time in our simulations the kinetic and potential energies are already of the same order of magnitude regardless the employed bottom motion (see Figures 9 & 10). Moreover, both curves $K(t)$ and $\Pi(t)$ continue to tend to the equilibrium state according to the theoretical predictions of the water wave theory, [96].

Despite some local attenuation effects of horizontal displacements on the free surface elevation, the complete generation scenario produces a tsunami wave with more important energy content. More precisely, our computations show that the horizontal displacements contribute about 10% into the total tsunami energy balance in the USGS scenario (this value is consistent with our previous results concerning the differences in wave amplitudes in Figure 8). This result is even more flagrant for the Caltech version which ascribes 24% of the energy to horizontal displacements. The free surface amplitudes in this case should differ as well by the same order of magnitude. As we already noted, the difference between trigonometric and exponential scenarios is negligible.
<figure><img src="image_17.png" /><figcaption>FIGURE 11. The total energy evolution predicted according to two different versions of the finite fault inversion. Note the different vertical scales on left (4.5 × 10¹²) and right (2.5 × 10¹³) images. The time t is given in seconds.</figcaption></figure>

Up to now the contribution of horizontal displacements has been quantified in terms of the bottom deformation (Figure 5), free surface elevation at a particular moment of time (t = 220 s, cf. Figure 8) and finally in terms of the wave energy (Figures 9 – 11). However, the wave energy cannot be directly measured in practice. Consequently, we continue to quantify the differences between complete and incomplete approaches in terms of the waveforms. Namely, we compute also the following relative measures of the waveforms difference which act on the whole computed free surface:

$$d_2(t) := \frac{||\eta_c(\vec{x}, t) - \eta_i(\vec{x}, t)||_{L_2(\mathbb{R}^2)}}{||\eta_i(\vec{x}, t)||_{L_2(\mathbb{R}^2)}}, \qquad d_{\infty}(t) := \frac{||\eta_c(\vec{x}, t) - \eta_i(\vec{x}, t)||_{L_{\infty}(\mathbb{R}^2)}}{||\eta_i(\vec{x}, t)||_{L_{\infty}(\mathbb{R}^2)}}. \qquad (3.1)$$

The simulation results are presented on Figure 12. One can see that both measures grow up to 15% and then oscillate around this level, at least during the simulation time. This result is in agreement with previous measurements of the horizontal displacements contribution based on the wave energy and the bottom deformation which were of the order of 10%. We note also that the curve $d_2(t)$ has a more regular behaviour than $d_{\infty}(t)$ since it is based on integral characteristics of the difference while the latter focusses on the characteristics of the extreme values.

It is also interesting to compare the tsunami energy with the energy of the underlying seismic event. The USGS Energy and Broadband solution indicates that the radiated seismic energy of the Java 2006 earthquake is equal to 3.2 × 10¹⁴ J. Hence, according to our computations with the USGS complete generation approach and the trigonometric scenario, about 0.68% of the seismic energy was transmitted to the tsunami wave (this portion raises to 7.14% for the Caltech version). The ratio of the total tsunami and seismic radiated energies can be used as a measure of the tsunami generation efficiency of a specific earthquake. The seismic radiated energy may not be the most appropriate
<figure><img src="image_18.png" /><figcaption>FIGURE 12. Relative difference between free surface elevations computed according to the complete and incomplete scenarios. For more details see equation (3.1). The horizontal axis is given in seconds, while the vertical scale is a percentage. The blue solid line (—) corresponds to the measure $d_2(t)$. The black dashed line (---) represents $d_\infty(t)$.</figcaption></figure>

parameter to use, but it has an advantage to be relatively robust and easily observable in contrast to the more relevant seismic fracture energy, [90, 51].

This parameter can also be estimated for the Great Sumatra-Andaman earthquake of December 26, 2004. The total energy of the Great Indian Ocean tsunami 2004 was estimated by T. Lay *et al*. (2005), [57], to be equal to $4.2 \times 10^{15}$ J. According to the same USGS Energy and Broadband solution the radiated seismic energy was equal to $1.1 \times 10^{17}$ J. Thus, the tsunamigenic efficiency of the Great Sumatra-Andaman earthquake is 3.81% which is bigger than that of Java 2006 event but remains in the same order of magnitude around 1%. It is possible that we do not still take into account all important factors in Java 2006 tsunami genesis. More examples of transmission of seismic energy to tsunami energy can be found e.g. in [58].

4. DISCUSSION AND CONCLUSIONS

In the present study the process of tsunami generation is further investigated. The current tsunamigenesis model relies on a combination of the finite fault solution, [45], and a recently proposed Weakly Nonlinear (WN) solver for the water wave problem with moving bottom, [24]. Consequently, in our model we incorporate recent advances in seismology and computational hydrodynamics.

This study is focused on the role of the horizontal displacements in the real world tsunami genesis process. By our intuition we know that a horizontal motion of the flat bottom will not cause any significant disturbance on the free surface. However, the real bathymetry
is far from being flat. The question which arises naturally is how to quantify the effect of
horizontal co-seismic displacements during real world events.

The primary goal of this study was to propose relatively simple, efficient and accurate procedures to model tsunami generation process in realistic environments. Special emphasis was payed to the role of horizontal displacements which should be also taken into account. Thus, the kinematics of horizontal co-seismic displacements were reconstructed and their effect on free surface motion was quantified. The evolution of kinetic and potential energies were also investigated in our study. In the case of July 17, 2006 Java event our simulations indicate that 10% of the energy input can be seemingly ascribed to effects of the horizontal bottom motion. This portion increases considerably if we switch to the scenario proposed by Caltech Tectonics Observatory [70].

The results presented in this study do not still explain the reasons for the extreme run-up values caused by July 17, 2006 Java tsunami, [26]. At least we hope that the proposed methodology illustrated on this important real world event will be proved helpful in future studies. In our opinion a successful theory should incorporate also other generation mechanisms such as local landslides/slumps which are subject to large uncertainties at the current stage of our understanding. Until now there are no detailed images of the seabed which could support or disprove this assumption. In this study we succeeded to quantify the significance of horizontal displacements for the tsunami generation.

APPENDIX A. APPLICATIONS TO SOME RECENT TSUNAMI EVENTS

In this Appendix we apply the techniques described in this manuscript to two recent significant events. We compute the energy transmission from the corresponding seismic event to the resulting tsunami wave with and without horizontal displacements. Throughout this Appendix we use the trigonometric scenario and the finite fault solutions produced by USGS.

**Remark 8.** In later versions of the finite fault solution the activation and rise time, as well as the seismic moment are specified for each sub-fault. Consequently, below we use this information to produce more accurate bottom kinematics. It allows us to compute also the seismic moment rate function<fnref n="5" />.

A.1. **Mentawai 2010 tsunami.** In October 25, 2010 a small portion of the subduction zone seaward of the Mentawai islands was ruptured by an earthquake of magnitude $M_w = 7.7$, [56]. This earthquake generated a tsunami wave with run-up values ranging from 3 to 9 meters (with even larger values at some places). On Pagai islands this tsunami caused more than 400 victims. This earthquake is characterized by 10° dip angle, a slow rupture velocity (≈ 1.5 km/s) which propagated during about 100 s over 100 km long source region. For our simulation we used the finite fault solution produced by USGS [38]. On Figure 13 we show the evolution of the total tsunami energy to be compared with the seismic moment release rate function.

<footnote n="5">⁵The total released seismic moment is given by the integral of this function over the rupture duration (time).</footnote>
<figure><img src="image_19.png" /><figcaption>FIGURE 13. Total tsunami energy and seismic moment rate function for the Mentawai 2010 event. The time t is given in seconds.</figcaption></figure>

Our computations give the following estimations:

* Complete scenario: $E_c = 1.11 \times 10^{12}$ J,

* Incomplete scenario: $E_i = 1.06 \times 10^{12}$ J.

Consequently, for this event only 4.5% of the total energy is due to horizontal displacements. According to the USGS energy and broadband solution the radiated seismic energy is $E_s = 1.4 \times 10^{15}$ J. The tsunami energy constitutes 0.08% of the radiated seismic energy. This surprising result can be explained by the fact that a big portion of co-seismic displacements occured on the land which did not fortunately contribute to the tsunami generation process.

### A.2. March 11, 2011 Japan earthquake and tsunami.
This tsunami event was caused by an undersea megathrust earthquake of magnitude $M_w = 9.0$ off the coast of Japan. Officially this earthquake was named the *Great East Japan Earthquake*. The Japanese National Police Agency has confirmed more than 13,000 deaths caused both by the earthquake and especially by the tsunami. Entire towns were devastated. The local infrastructure including the Fukushima Nuclear Power Plant were heavily affected with consequences which are widely known.

We use again the finite fault inversion performed at the USGS [39]. The tsunami total energy evolution along with the moment rate function are represented on Figure 14. At the end of the simulation we obtain the following result:

* Complete scenario: $E_c = 1.64 \times 10^{15}$ J,

* Incomplete scenario: $E_i = 1.50 \times 10^{15}$ J.

Henceforth, the contribution of horizontal displacements is estimated to be 9.2%. The radiated seismic energy is $E_s = 1.9 \times 10^{17}$ J according to the USGS energy and broadband solution. The total tsunami wave energy represents only 0.87% of the radiated seismic
<figure><img src="image_20.png" /><figcaption>FIGURE 14. Total tsunami energy and seismic moment rate function for the Great East Japan Earthquake. The time *t* is given in seconds.</figcaption></figure>

energy. This value is much lower than the corresponding value for the Great Sumatra-Andaman earthquake of December 26, 2004.

## ACKNOWLEDGEMENT

D. Dutykh acknowledges the support from French Agence Nationale de la Recherche, project MathOcean (Grant ANR-08-BLAN-0301-01) and Ulysses Program of the French Ministry of Foreign Affairs under the project 23725ZA. French-Russian collaboration is supported by CNRS PICS project No. 5607. L. Chubarov and Yu. Shokin acknowledge the support from the Russian Foundation for Basic Research under RFBR Project No. 10-05-91052 and from the Basic Program of Fundamental Research SB RAS Project No. IV.31.2.1.

We would like to thank Professors Didier Clamond and Frédéric Dias for very helpful discussions on numerical simulation of water waves. Special thanks go to Professor Costas Synolakis whose work on tsunami waves has also been the source of our inspiration. Finally we thank two anonymous referees for the time they spent to review our manuscript and whose comments helped us a lot to improve its quality.

## REFERENCES

[1] C. J. Ammon, C. Ji, H.-K. Thio, D. Robinson, S. Ni, V. Hjorleifsdottir, H. Kanamori, T. Lay, S. Das, D. Helmberger, G. Ichinose, J. Polet, and D. Wald. Rupture process of the 2004 Sumatra-Andaman earthquake. *Science*, 308:1133-1139, 2005. 1

[2] C. J. Ammon, H. Kanamori, T. Lay, and A. A. Velasco. The 17 July 2006 Java tsunami earthquake. *Geophysical Research Letters*, 33:L24308, 2006. 2, 13

[3] T. Baba, P. R. Cummins, T. Hori, and Y. Kaneda. High precision slip distribution of the 1944 Tonankai earthquake inferred from tsunami waveforms: Possible slip on a splay fault. *Tectonophysics*, 426:119-134, 2006. 9
[4] J.-P. Bardet, C. E. Synolakis, H. L. Davies, F. Imamura, and E. A. Okal. Landslide Tsunamis: Recent Findings and Research Directions. *Pure Appl. Geophys.*, 160:1793–1809, 2003. 3

[5] C. Bassin, G. Laske, and G. Masters. The Current Limits of Resolution for Surface Wave Tomography in North America. *EOS Trans AGU*, 81:F897, 2000. 5

[6] S. Beisel, L. Chubarov, I. Didenkulova, E. Kit, A. Levin, E. Pelinovsky, Yu. Shokin, and M. Sladkevich. The 1956 Greek tsunami recorded at Yafo, Israel, and its numerical modeling. *Journal of Geophysical Research - Oceans*, 114:C09002, 2009. 9

[7] E. Berg. Field survey of the Tsunamis of 28 March 1964 in Alaska, and conclusions as to the origin of the major Tsunami. Technical report, Hawaii Institute of Geophysics, University of Hawaii, Honolulu, 1970. 3, 5

[8] R. Bilham. A flying start, then a slow slip. *Science*, 308:1126–1127, 2005. 2

[9] D. Clamond and D. Dutykh. Practical use of variational principles for modeling water waves. *Physica D: Nonlinear Phenomena*, 241(1):25–36, 2012. 13

[10] R. R. Coifman and Y. Meyer. Nonlinear harmonic analysis and analytic dependence. *Pseudodifferential Operators and Applications*, 43:71–78, 1985. 12

[11] W. Craig and C. Sulem. Numerical simulation of gravity waves. *J. Comput. Phys.*, 108:73–83, 1993. 11, 12

[12] W. Craig, C. Sulem, and P.-L. Sulem. Nonlinear modulation of gravity waves: a rigorous approach. *Nonlinearity*, 5(2):497–522, 1992. 11

[13] S. Das and P. Suhadolc. On the inverse problem for earthquake rupture: The Haskell-type source model. *Journal of Geophysical Research*, 101(B3):5725–5738, 1996. 8

[14] A. I. Delis, M. Kazolea, and N. A. Kampanis. A robust high-resolution finite volume scheme for the simulation of long waves over complex domains. *Int. J. Numer. Meth. Fluids*, 56:419–452, 2008. 4

[15] F. Dias and T. J. Bridges. The numerical computation of freely propagating time-dependent irrotational water waves. *Fluid Dynamics Research*, 38:803–830, 2006. 11

[16] F. Dias and D. Dutykh. *Dynamics of tsunami waves*, chapter Dynamics of pages 35–60. Springer, 2007. 4, 7, 10

[17] J. R. Dormand and P. J. Prince. A family of embedded Runge-Kutta formulae. *J. Comp. Appl. Math.*, 6:19–26, 1980. 14

[18] D. Dutykh. *Mathematical modelling of tsunami waves*. Ph.d. thesis, École Normale Supérieure de Cachan, December 2007. 3, 7, 10

[19] D. Dutykh and F. Dias. Water waves generated by a moving bottom. In Anjan Kundu, editor, *Tsunami and Nonlinear waves*, pages 65–96. Springer Verlag (Geo Sc.), 2007. 3, 5, 7, 8, 10

[20] D. Dutykh and F. Dias. Energy of tsunami waves generated by bottom motion. *Proc. R. Soc. A*, 465:725–744, 2009. 4, 12

[21] D. Dutykh and F. Dias. Tsunami generation by dynamic displacement of sea bed due to dip-slip faulting. *Mathematics and Computers in Simulation*, 80(4):837–848, 2009. 3, 5

[22] D. Dutykh, F. Dias, and Y. Kervella. Linear theory of wave generation by a moving bottom. *C. R. Acad. Sci. Paris, Ser. I*, 343:499–504, 2006. 8, 10, 16

[23] D. Dutykh, Th. Katsaounis, and D. Mitsotakis. Finite volume schemes for dispersive wave propagation and runup. *J. Comput. Phys*, 230:3035–3061, 2011. 4

[24] D. Dutykh, D. Mitsotakis, X. Gardeil, and F. Dias. On the use of the finite fault solution for tsunami generation problems. *Theor. Comput. Fluid Dyn.*, March 2012. 4, 5, 6, 7, 12, 16, 22

[25] D. Dutykh, R. Poncet, and F. Dias. The VOLNA code for the numerical modeling of tsunami waves: Generation, propagation and inundation. *Eur. J. Mech. B/Fluids*, 30(6):598–615, 2011. 3, 4

[26] H. M. Fritz, W. Kongko, A. Moore, B. McAdoo, J. Goff, C. Harbitz, B. Uslu, N. Kalligeris, D. Suteja, K. Kalsum, V. V. Titov, A. Gusman, H. Latief, E. Santoso, S. Sujoko, D. Djulkarnaen, H. Sunendar, and C. Synolakis. Extreme runup from the 17 July 2006 Java tsunami. *Geophys. Res. Lett.*, 34:L12602, 2007. 2, 13, 23
[27] D. Fructus, D. Clamond, O. Kristiansen, and J. Grue. An efficient model for threedimensional surface wave simulations. Part I: Free space problems. *J. Comput. Phys.*, 205:665–685, 2005. 11

[28] D. Fructus and J. Grue. An explicit method for the nonlinear interaction between water waves and variable and moving bottom topography. *J. Comp. Phys.*, 222:720–739, 2007. 5

[29] D. R. Fuhrman and P. A. Madsen. Tsunami generation, propagation, and run-up with a high-order Boussinesq model. *Coastal Engineering*, 56(7):747–758, 2009. 5

[30] Y. Fujii and K. Satake. Source of the July 2006 West Java tsunami estimated from tide gauge records. *Geophys. Res. Lett.*, 33:L24317, 2006. 6

[31] E. Geist. Complex earthquake rupture and local tsunamis. *J. Geophys. Res.*, 107:B5, 2002. 5

[32] E. L. Geist, S. L. Bilek, D. Arcas, and V. V. Titov. Differences in tsunami generation between the December 26, 2004 and March 28, 2005 Sumatra earthquakes. *Earth Planets Space*, 58:185–193, 2006. 2, 5

[33] V. K. Gusiakov. Residual displacements on the surface of elastic half-space. In *Conditionally correct problems of mathematical physics in interpretation of geophysical observations*, pages 23–51, 1978. 3

[34] P. Guyenne and D. P. Nicholls. A high-order spectral method for nonlinear water waves over moving bottom topography. *SIAM J. Sci. Comput.*, 30(1):81–101, 2007. 12

[35] J. Hammack. A note on tsunamis: their generation and propagation in an ocean of uniform depth. *J. Fluid Mech.*, 60:769–799, 1973. 8, 9, 16

[36] J. L. Hammack. *Tsunamis - A Model of Their Generation and Propagation*. PhD thesis, California Institute of Technology, 1972. 8

[37] C. Harbitz. Model simulations of tsunamis generated by the Storegga Slides. *Marine Geology*, 104:1–21, 1992. 3

[38] G. Hayes. Finite Fault Model Result of the Oct 25, 2010 Mw 7.7 Southern Sumatra Earthquake. Technical report, USGS, 2010. 23

[39] G. Hayes. Finite Fault Model Result of the Mar 11, 2011 Mw 9.0 Earthquake Offshore Honshu, Japan. Technical report, USGS, 2011. 24

[40] G. Ichinose, P. Somerville, H. K. Thio, R. Graves, and D. O'Connell. Rupture process of the 1964 Prince William Sound, Alaska, earthquake from the combined inversion of seismic, tsunami, and geodetic data. *Journal of Geophysical Research*, 112(B7):B07306, July 2007. 6

[41] G. A. Ichinose, K. Satake, J. G. Anderson, and M. M. Lahren. The potential hazard from tsunami and seiche waves generated by large earthquakes within Lake Tahoe, California-Nevada. *Geophys. Res. Lett.*, 27:1203–1206, 2000. 3

[42] F. Imamura. *Long-wave runup models*, chapter Simulation, pages 231–241. World Scientific, 1996. 4, 5

[43] F. Imamura, A. C. Yalciner, and G. Ozyurt. *Tsunami modelling manual*, April 2006. 5

[44] C. Ji. Preliminary Result of the 2006 July 17 Magnitude 7.7 - South of Java, Indonesia Earthquake. Technical report, USGS, 2006. 6, 7, 8, 13

[45] C. Ji, D. J. Wald, and D. V. Helmberger. Source description of the 1999 Hector Mine, California earthquake; Part I: Wavelet domain inversion theory and resolution analysis. *Bull. Seism. Soc. Am.*, 92(4):1192–1207, 2002. 5, 22

[46] J. M. Johnson and K. Satake. Source parameters of the 1957 Aleutian earthquake from tsunami waveforms. *Geophysical Research Letters*, 20:1487–1490, 1993. 2

[47] K. Kajiura. The leading wave of tsunami. *Bull. Earthquake Res. Inst., Tokyo Univ.*, 41:535–571, 1963. 2, 4

[48] K. Kajiura. Tsunami source, energy and the directivity of wave radiation. *Bull. Earthquake Research Institute*, 48:835–869, 1970. 5

[49] H. Kanamori. The Alaska earthquake of 1964: Radiation of Long-Period Surface Waves and Source Mechanism. *Journal of Geophysical Research*, 75(26):5029–5040, 1970. 2
[50] H. Kanamori. Mechanism of tsunami earthquakes. *Physics of the Earth and Planetary Interiors*, 6(5):346–359, January 1972. 8

[51] H. Kanamori and L. Rivera. Energy partitioning during an earthquake. *Geophysical Monograph Series*, 170:3–13, 2006. 22

[52] Y. Kervella, D. Dutykh, and F. Dias. Comparison between three-dimensional linear and nonlinear tsunami generation models. *Theor. Comput. Fluid Dyn.*, 21:245–269, 2007. 5, 8, 10, 12, 16

[53] D.-H. Kim, Y.-S. Cho, and Y.-K. Yi. Propagation and run-up of nearshore tsunamis with HLLC approximate Riemann solver. *Ocean Engineering*, 34:1164–1173, 2007. 4

[54] H. Lamb. *Hydrodynamics*. Cambridge University Press, 1932. 10

[55] J. F. Lander. Tsunamis affecting Alaska 1737 - 1996. *U.S. Department of Commerce, National Oceanic and Atmospheric Administration*, 31:195, 1996. 2

[56] T. Lay, C. J. Ammon, H. Kanamori, Y. Yamazaki, K. F. Cheung, and A. R. Hutko. The 25 October 2010 Mentawai tsunami earthquake (Mw 7.8) and the tsunami hazard presented by shallow megathrust ruptures. *Geophys. Res. Lett.*, 38:L06302, 2011. 23

[57] T. Lay, H. Kanamori, C. J. Ammon, M. Nettles, S. N. Ward, R. C. Aster, S. L. Beck, S. L. Bilek, M. R. Brudzinski, R. Butler, H. R. DeShon, G. Ekstrom, K. Satake, and S. Sipkin. The great Sumatra-Andaman earthquake of 26 December 2004. *Science*, 308:1127–1133, 2005. 2, 22

[58] F. Løvholt, H. Bungum, C. B. Harbitz, S. Glimsdal, C. Lindholm, and G. Pedersen. Earthquake related tsunami hazard along the western coast of Thailand. *Nat. Hazards Earth Syst. Sci.*, 6:1–18, 2006. 22

[59] L. Mansinha and D. E. Smylie. Effect of earthquakes on the Chandler wobble and the secular polar shift. *J. Geophys. Res.*, 72:4731–4743, 1967. 3, 7

[60] L. Mansinha and D. E. Smylie. The displacement fields of inclined faults. *Bull. Seism. Soc. Am.*, 61:1433–1440, 1971. 3, 7

[61] J. McCloskey, A. Antonioli, A. Piatanesi, K. Sieh, S. Steacy, M. Nalbant S. Cocco, C. Giunchi, J. D. Huang, and P. Dunlop. Tsunami threat in the Indian Ocean from a future megathrust earthquake west of Sumatra. *Earth and Planetary Science Letters*, 265:61–81, 2008. 5

[62] C. C. Mei. *The applied dynamics of ocean surface waves*. World Scientific, 1994. 10

[63] R. D. Mindlin. Force at a point in the interior of a semi-infinite medium. *Physics*, 7:195–202, 1936. 7

[64] D. E. Mitsotakis. Boussinesq systems in two space dimensions over a variable bottom for the generation and propagation of tsunami waves. *Math. Comp. Simul.*, 80:860–873, 2009. 3, 16

[65] T. Ohmachi, H. Tsukiyama, and H. Matsumoto. Simulation of tsunami induced by dynamic displacement of seabed due to seismic faulting. *Bull. Seism. Soc. Am.*, 91:1898–1909, 2001. 3

[66] Y. Okada. Surface deformation due to shear and tensile faults in a half-space. *Bull. Seism. Soc. Am.*, 75:1135–1154, 1985. 2, 4, 7

[67] Y. Okada. Internal deformation due to shear and tensile faults in a half-space. *Bull. Seism. Soc. Am.*, 82:1018–1040, 1992. 2, 4, 7

[68] E. Okal and C. E. Synolakis. Field survey and numerical simulations: a theoretical comparison of tsunamis from dislocations and landslides. *Pure Appl. Geophys.*, 160:2177–2188, 2003. 2

[69] E. A. Okal and C. E. Synolakis. A theoretical comparison of tsunamis from dislocations and landslides. *Pure and Applied Geophysics*, 160:2177–2188, 2003. 3

[70] A. Ozgun Konca. Preliminary Result 06/07/17 (Mw 7.9), Southern Java Earthquake. Technical report, http://www.tectonics.caltech.edu/slip_history/2006, USGS, 2006. 6, 13, 19, 20, 23

[71] A. Piatanesi, S. Tinti, and G. Pagnoni. Tsunami waveform inversion by numerical finite-elements Green's functions. *Nat. Hazards Earth Syst. Sci.*, 1:187–194, 2001. 6

[72] C. Pires and P. M. A. Miranda. Sensitivity of the adjoint method in the inversion of tsunami source parameters. *Nat. Hazards Earth Syst. Sci.*, 3:341–351, 2003. 6

[73] F. Press. Displacements, strains and tilts at tele-seismic distances. *J. Geophys. Res.*, 70:2395–2412, 1965. 7
[74] A. B. Rabinovich, L. I. Lobkovsky, I. V. Fine, R. E. Thomson, T. N. Ivelskaya, and E. A. Kulikov. Near-source observations and modeling of the Kuril Islands tsunamis of 15 November 2006 and 13 January 2007. *Adv. Geosci.*, 14:105–116, 2008. 3

[75] J. Rhie, D. Dreger, R. Burgmann, and B. Romanowicz. Slip of the 2004 Sumatra-Andaman Earthquake from Joint Inversion of Long-Period Global Seismic Waveforms and GPS Static Offsets. *Bulletin of the Seismological Society of America*, 97(1A):S115-S127, January 2007. 6

[76] T. Saito and T. Furumura. Three-dimensional tsunami generation simulation due to sea-bottom deformation and its interpretation based on the linear theory. *Geophys. J. Int.*, 178:877–888, 2009. 3, 12, 16

[77] K. Satake. Linear and nonlinear computations of the 1992 Nicaragua earthquake tsunami. *Pure Appl. Geophys.*, 144(3-4):455–470, 1995. 12

[78] Y. T. Song, L.-L. Fu, V. Zlotnicki, C. Ji, V. Hjorleifsdottir, C. K. Shum, and Y. Yi. The role of horizontal impulses of the faulting continental slope in generating the 26 December 2004 tsunami. *Ocean Modelling*, 20:362–379, 2008. 3, 4, 5, 9

[79] Y. T. Song and T. Y. Hou. Parametric vertical coordinate formulation for multiscale, Boussinesq, and non-Boussinesq ocean modeling. *Ocean Modelling*, 11(3-4):298–332, January 2006. 12

[80] J. A. Steketee. On Volterra's dislocation in a semi-infinite elastic medium. *Can. J. Phys.*, 36:192–205, 1958. 4, 7

[81] D. J. Stevenson. Tsunamis and Earthquakes: What Physics is interesting? *Physics Today*, June:10–11, 2005. 3

[82] J. J. Stoker. *Water waves, the mathematical theory with applications*. Wiley, 1958. 10

[83] C. E. Synolakis and E. N. Bernard. Tsunami science before and beyond Boxing Day 2004. *Phil. Trans. R. Soc. A*, 364:2231–2265, 2006. 2, 4

[84] Y. Tanioka and K. Satake. Tsunami generation by horizontal displacement of ocean bottom. *Geo-physical Research Letters*, 23:861–864, 1996. 3, 4, 5, 9, 12

[85] V. V. Titov and F. I. González. Implementation and testing of the method of splitting tsunami (MOST) model. Technical Report ERL PMEL-112, Pacific Marine Environmental Laboratory, NOAA, 1997. 5

[86] M. I. Todorovska, A. Hayir, and M. D. Trifunac. A note on tsunami amplitudes above submarine slides and slumps. *Soil Dynamics and Earthquake Engineering*, 22:129–141, 2002. 8

[87] M. I. Todorovska and M. D. Trifunac. Generation of tsunamis by a slowly spreading uplift of the seafloor. *Soil Dynamics and Earthquake Engineering*, 21:151–167, 2001. 5, 8

[88] Y. Tsuji, F. Imamura, H. Matsutomi, C. E. Synolakis, P. T. Nanang, Jumadi, S Harada, S. S. Han, K. Arai, and B. Cook. Field survey of the East Java earthquake and tsunami of June 3, 1994. *Pure Appl. Geophys.*, 144:839–854, 1995. 2

[89] A. S. Velichko, S. F. Dotsenko, and E. N. Potetyunko. Amplitude-energy characteristics of tsunami waves for various types of seismic sources generating them. *Phys. Oceanogr.*, 12(6):308–322, 2002. 4

[90] A. Venkataraman and H. Kanamori. Observational constraints on the fracture energy of subduction zone earthquakes. *J. Geophys. Res.*, 109:B05302, 2004. 22

[91] J. K. Vennard and R. L. Street. *Elementary fluid mechanics*. Springer, sixth edition, 1982. 3, 4

[92] V. Volterra. Sur l'équilibre des corps élastiques multiplement connexes. *Annales Scientifiques de l'Ecole Normale Supérieure*, 24(3):401–517, 1907. 7

[93] K. T. Walker, M. Ishii, and P. M. Shearer. Rupture details of the 28 March 2005 Sumatra Mw 8.6 earthquake imaged with teleseismic P waves. *Geophysical Research Letters*, 32:24303, 2005. 2

[94] X. Wang and P.L.-F. Liu. An analysis of 2004 Sumatra earthquake fault plane mechanisms and Indian Ocean tsunami. *J. Hydr. Res.*, 44(2):147–154, 2006. 5

[95] S. N. Ward. Landslide tsunami. *J. Geophysical Res.*, 106:11201–11215, 2001. 3

[96] G. B. Whitham. *Linear and nonlinear waves*. John Wiley & Sons Inc., New York, 1999. 10, 20
[97] Y. Yagi. Source rupture process of the 2003 Tokachi-oki earthquake determined by joint inversion of teleseismic body wave and strong ground motion data. *Earth Planets Space*, 56(3):311–316, 2004. 6

[98] V. E. Zakharov. Stability of periodic waves of finite amplitude on the surface of a deep fluid. *J. Appl. Mech. Tech. Phys.*, 9:190–194, 1968. 11

LAMA, UMR 5127 CNRS, UNIVERSITÉ DE SAVOIE, CAMPUS SCIENTIFIQUE, 73376 LE BOURGET-DU-LAC CEDEX, FRANCE

*E-mail address:* Denys.Dutykh@univ-savoie.fr

URL: http://www.lama.univ-savoie.fr/~dutykh/

IMA, UNIVERSITY OF MINNESOTA, 114 LIND HALL, 207 CHURCH STREET SE, MINNEAPOLIS MN 55455, USA

*E-mail address:* dmitsot@gmail.com

URL: http://sites.google.com/site/dmitsot/

INSTITUTE OF COMPUTATIONAL TECHNOLOGIES, SIBERIAN BRANCH OF THE RUSSIAN ACADEMY OF SCIENCES, 6 ACAD. LAVRENTJEV AVENUE, 630090 NOVOSIBIRSK, RUSSIA

*E-mail address:* chubarov@ict.nsc.ru

URL: http://www.ict.nsc.ru/

INSTITUTE OF COMPUTATIONAL TECHNOLOGIES, SIBERIAN BRANCH OF THE RUSSIAN ACADEMY OF SCIENCES, 6 ACAD. LAVRENTJEV AVENUE, 630090 NOVOSIBIRSK, RUSSIA

*E-mail address:* shokin@ict.nsc.ru

URL: http://www.ict.nsc.ru/