## 1. Introduction

The coastal areas of Tohoku district suffered serious damage from the tsunami caused by the 2011 off the Pacific Coast of Tohoku Earthquake that occurred on March 11, 2011 [1, 2]. Numerical simulations are used to develop disaster-prevention measures to deal with such tsunami disasters. They are also used to predict potential future tsunami disasters, to design disaster-prevention facilities such as coastal breakwaters and levees and to predict tsunami attacks immediately after an earthquake occurs [3, 4]. The Central Disaster Prevention Council [5] prepares a basic disaster prevention plan and participates in the determination of important disaster-prevention matters. If the Tonankai-Nankai earthquake occurs, it may cause considerable damage. Consequently, the council has been performing numerical calculations to predict wave height and arrival time when tsunami reaches the coast.

In this chapter, to let it be self-contained, the viscous shallow-water equations are again derived from the Navier-Stokes equations in which the hydrostatic pressure in the direction of gravity is assumed. In the numerical analysis of tsunamis, the viscosity term is often omitted or simply added [6, 7]. In this study, however, a computational method that does not omit the viscosity term is adopted, that is, because it enables more rigorous analysis to be performed and we intend to include the viscosity term in future stress analysis of tsunamis. The approximation scheme [8] is given below and simulation results [9] for a Tohoku-Oki model are again presented for self-containing. Our main concern in this chapter is how to set the boundary condition on the open boundary, which is completely changed to a new one.

## 2. Viscous shallow-water equations

The tsunami generated by the GEJE caused serious damage to the coastal areas of the Tohoku district. Numerical simulations are used to predict damage caused by tsunamis. Shallow-water equations are generally used in numerical simulations of tsunami propagation from the open sea to the coast. This subsection focuses on viscous shallow-water equations and attempts to generate a computational method using finite-element techniques based on the previous investigations of Kanayama and Ohtsuka [8].

In the numerical analysis of tsunami, a viscosity term is often omitted or simply added [6, 7]. In this subsection, however, a computational method that does not omit the viscosity term is adopted, that is, because it enables more rigorous analysis to be performed, and the authors intend to include the viscosity term in future stress analysis of tsunami. Recently, the authors have found an interesting book [10] that deals with sound mathematical topics related to the viscous shallow-water equations.

### 2.1. Derivation of the viscous shallow-water equations

As shown in **Figure 1**, Kanayama and Dan [9] considered the shallow-water long-wave flow in which the wavelength is sufficiently long relative to the water depth [11]. First, the viscous shallow-water equations are again derived to let this chapter be self-contained. Detailed derivation comes from the study of Kanayama and Ushijima [11]. Orthogonal coordinates [m] are used, where *x*_{1} and *x*_{2} represent directions in the horizontal plane and *x*_{3} represents the vertical direction, and the time [s] is represented by *t*. Assuming the hydrostatic pressure in the *x*_{3} direction, the following Navier-Stokes equations are used with the external forces, namely the Coriolis forces and the gravity are assumed to act in the *x*_{1} and *x*_{2} directions and the *x*_{3} direction, respectively.

*u*(*x*_{1}, *x*_{2}, *x*_{3}, *t*) represents the fluid velocity component [m/s] in the *x _{i}* (

*i*=1−3) direction,

*p*(

*x*

_{1},

*x*

_{2},

*x*

_{3},

*t*) denotes the pressure [N/m

^{2}],

*ρ*is the density [kg/m

^{3}],

*τ*is the stress component [N/m

_{ij}^{2}] in the

*x*direction acting on the

_{i}*x*plane,

_{j}*f*is the Coriolis coefficient [1/s], and

*g*is the acceleration [m/s

^{2}] due to the gravity. In addition, the ordinate [m] of the water surface (water level) is represented by

*ζ*and the ordinate [m] of the bottom surface is represented by

*h*. Note that

*h*usually has negative input values depending on

*x*(

_{i}*i*= 1–2).

Eqs. (1)–(3) [9] are integrated with respect to *x*_{3} from the bottom *h* to the water surface *ζ*. The velocity component in the normal direction is assumed to be zero at the water and bottom surfaces. By doing this, the viscous shallow-water equations expressed by the water level and the averaged velocity can be derived. If the layer thickness is *H*(*x*_{1}, *x*_{2}, *t*) and the average velocity in the *x _{i}* (

*i*= 1, 2) direction is

*U*(

_{i}*x*

_{1},

*x*

_{2},

*t*), then

*H*and

*U*are given by the following equations:

_{i}Note that *H* is assumed to be positive in this chapter, whose guarantee is a difficult mathematical problem [9].

The fixed density of the layer is represented by *ρ* and the stress component in the *x*_{i} direction acting on the *x _{j}* plane is represented by

^{2}] is represented by

*μ*

_{H}, the wind effect coefficient is represented by

*θ*, the air density is represented by

*ρ*, the wind velocity component speed in the

_{a}*x*direction is represented by

_{i}*W*, and the Chezy coefficient [m

_{i}^{1/2}/s] is represented by

*C*. Then, the following viscous shallow-water equations are derived [9]:

(7) |

where

Note that

Finite-element approximation is performed for the viscous shallow-water Eqs. (6) and (7) with given initial conditions and boundary conditions.

The computational domain is the two-dimensional (2D) polygonal region *Ω* surrounded by the boundaries *Γ*_{c} and *Γ*_{o}. In this region, the orthogonal coordinates *x* = (*x*_{1}, *x*_{2}) are used. The boundary conditions on *Γ*_{c }and *Γ*_{o} and the initial conditions at *t* = 0 are as follows [9]:

Boundary conditions, Eqs. (9) and (10):

Here, *n _{i}* is the component in the

*x*direction of the unit outward normal vector on the boundary. On

_{i}*Γ*

_{o},

Initial conditions:

Here, *U*_{i0}(*x*) and *ζ*_{0}(*x*) are the initial values of *U*_{i} and *ζ*, respectively.

### 2.2. Finite-element approximations for the viscous shallow-water equations

The finite-element approximation [8] is as follows. First, Eqs. (6) and (7) are multiplied by test functions and then integrated over the computational domain *Ω*. Subsequently, the approximation is carried out for terms containing a derivative with respect to time by using an explicit method. For terms containing spatial derivatives, the finite-element approximation is carried out by using the piecewise linear basis function *k* and is 0 at other places.

(13) |

In the above, *k* after *n* time steps, and Δ*t* represents the size of time steps. In addition, the following symbols are also used:

It is well known that integration by parts is used for the viscosity term in (13).

Because of the presence of the basic boundary conditions, Eq. (12) holds for all nodes except for the nodes on *Γ*_{o} and Eq. (13) holds for all the nodes except for the nodes on *Γ*_{c}. However, on the boundary *Γ*_{o}, an approximation method is adopted in which the advective term of Eq. (13) uses Tabata’s upwind approximation [9, 12]. A mathematical justification for the approximation scheme for linearized equations related to the above scheme was given by Kanayama and Ushijima [13, 14].

In general, tsunami is excited in the following two ways. The first one is to consider the tsunami excitation as the initial condition of the water surface, for which we do not have sufficient input information in such an artificial tsunami of Hakata Bay [9]. The second one is to consider it as the boundary condition of the water surface as in the next subsection. In the setting of Hakata Bay, a computational domain is not so wide that the above approach may be the only way. It is also noted that 50 [m] at the open boundary for the later Tohoku-Oki case may be too high. In the computation, the tsunami arrived at Oshika Peninsula after 20 min, and the highest wave height reached 15 [m]. These numerical results should be checked more carefully with data on the open boundary.

## 3. Computational examples

### 3.1. Tohoku-Oki models

Now, numerical results [9] for a Tohoku-Oki model are again shown for self-containing. When not specifically defined, the same physical values as for Hakata Bay were used. Regarding the mesh data, the total number of nodes is 56,562, the total number of elements is 113,013, and the ocean floor ordinate is set to become deep gradually from *h* = ‒10 to *h* = ‒1,000 [m]. We set the initial conditions at ζ_{0} (*x*) = 0 and *U*_{i0} (*x*) = 0, and the boundary conditions *t* < 60 [s]) at the tsunami-generating area and **Figure 2**). It is noted that 50 [m] at *Γ*_{o} may be too high. In our computation, the tsunami arrived at Oshika Peninsula in **Figure 2** after 20 min and the highest wave height reached 15 [m]. These numerical results should more carefully be checked with boundary data on *Γ*_{o}.

**Figure 2a** shows the ordinates of the water surface (water level) after 18 [min]. Since the computational domain is not wide, it looks that there is a reflection from the northern boundary in **Figure 2a**. This artificial reflection can be removed by suitable boundary conditions on *Γ*_{o}. Details are later mentioned. **Figure 2b** shows the water level change at the two points A and B in **Figure 2**. The wave height at the point B is higher than the point A after about 1500 [s]. When the tsunami reaches coastal areas, the water depth becomes shallow and the wave height becomes high.

Next, new numerical results for the second Tohoku-Oki model are shown. When not specifically defined, the same physical values as for Hakata Bay were used. Regarding the mesh data, the total number of nodes is 563,100, the total number of elements is 1,123,178, and **Figure 3** shows the ocean floor ordinate, which is constructed from 30-arc sec interval grid of JTOPO30, provided by the Marine Information Research Center. We set the initial conditions at ζ_{0} (*x*) like **Figure 4** based on the data of Fujii et al. [15] as a source of tsunami. **Figure 5a–f** shows the ordinates of the water surface (water level) every 8 [min] from after 8 [min] to after 48 [min]. In our previous paper [9], since the computational domain is not wide, it looks that there is a reflection from the boundary *Γ*_{o} (see **Figure 2a**). This artificial reflection can be removed by changing boundary conditions from (10) to (16) as follows:

where *U*_{n} and *U*_{t} denote the normal component of velocity and the tangential component of velocity, respectively, and *c* is a constant. We used (*c* = 0.9) in this chapter.

**Figure 6** shows the water level change at the three points Miyako, Soma and Choshi in **Figure 4**. Solid curves show observed data [15] and dashed curves show numerical results. Though the arrival time of the first wave is almost same, the wave height is small compared with observed data at the points of Miyako and Soma. Results may be improved by changing initial and boundary conditions.

## 4. Concluding remarks

In this study, the viscous shallow-water equations have again been derived from the Navier-Stokes equations in which the hydrostatic pressure in the direction of gravity is assumed. Tsunami propagation is then simulated by a finite-element computation. In the study of Kanayama and Dan [9] using the Hakata Bay model, in which the tsunami-generating area was taken to be the epicenter of the West off Fukuoka Prefecture Earthquake (2005), tsunami propagation from the open sea to the coastal area might be produced. This is an example in which the conventional open boundary condition Eq. (10) is still valid. In this chapter, however, in the second Tohoku-Oki model, the open boundary condition is changed from Eq. (10) to Eq. (16) to overcome artificial reflection. In addition, we have reconstructed a system in which the tsunami arrival time and the wave height at the time of tsunami attack can be obtained from numerical results. Since our analysis takes the viscosity term into account, this study can also be considered to be a preliminary study for future planned stress analyses of tsunamis (see **Figure 7**). We still want to check whether the hydrostatic pressure is always the dominant force in the case of tsunami attack compared with other forces including the viscous force.