**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.

Tsunami Propagation from the Open Sea to the Coast http://dx.doi.org/10.5772/63814 63

$$\sum\_{j=1}^{3} \frac{\partial u\_j}{\partial \mathbf{x}\_j} = \mathbf{0},\tag{1}$$

$$\frac{\partial \widehat{u}\_{1}}{\partial t} + \sum\_{j=1}^{3} u\_{j} \frac{\partial \widehat{u}\_{1}}{\partial \mathbf{x}\_{j}} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}\_{1}} + \frac{1}{\rho} \sum\_{j=1}^{3} \frac{\partial \pi\_{1j}}{\partial \mathbf{x}\_{j}} + f \widehat{u}\_{2} \,, \tag{2}$$

$$\frac{\partial u\_2}{\partial t} + \sum\_{j=1}^3 u\_j \frac{\partial u\_2}{\partial \mathbf{x}\_j} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}\_2} + \frac{1}{\rho} \sum\_{j=1}^3 \frac{\partial \mathbf{r}\_{2j}}{\partial \mathbf{x}\_j} - f \mathbf{u}\_1,\tag{3}$$

$$0 = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}\_3} - \mathbf{g}.\tag{4}$$

*u*(*x*1, *x*2, *x*3, *t*) represents the fluid velocity component [m/s] in the *xi* (*i*=1−3) direction, *p*(*x*1, *x*2, *x*3, *t*) denotes the pressure [N/m2 ], *ρ* is the density [kg/m3 ], *τij* is the stress component [N/m2 ] in the *xi* direction acting on the *xj* plane, *f* is the Coriolis coefficient [1/s], and *g* is the acceleration [m/s2 ] 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 *xi* (*i* = 1–2).

**Figure 1.** Variables of the computational model [9].

disaster-prevention matters. If the Tonankai-Nankai earthquake occurs, it may cause considerable damage. Consequently, the council has been performing numerical calculations to predict

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

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

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

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

condition on the open boundary, which is completely changed to a new one.

wave height and arrival time when tsunami reaches the coast.

62 Tsunami

**2. Viscous shallow-water equations**

investigations of Kanayama and Ohtsuka [8].

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

viscous shallow-water equations.

the *x*3 direction, respectively.

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 *xi* (*i* = 1, 2) direction is *Ui* (*x*1, *x*2, *t*), then *H* and *Ui* are given by the following equations:

$$H = \zeta - h, \ U\_i = \frac{1}{H} \int\_h^\zeta \mu\_i \, d\mathbf{x}\_3. \tag{5}$$

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*<sup>i</sup> direction acting on the *xj* plane is represented by . The horizontal viscosity constant [Ns/m2 ] is represented by *μ*H, the wind effect coefficient is represented by *θ*, the air density is represented by *ρa*, the wind velocity component speed in the *xi* direction is represented by *Wi* , and the Chezy coefficient [m1/2/s] is represented by *C*. Then, the following viscous shallow-water equations are derived [9]:

$$\frac{\partial \mathcal{L}}{\partial t} + \sum\_{j=1}^{2} \frac{\partial}{\partial \mathbf{x}\_{j}} \left( H U\_{j} \right) = 0,\tag{6}$$

$$\begin{split} &H\left(\frac{\partial U\_{i}}{\partial t} + \sum\_{j=1}^{2} U\_{j}\frac{\partial U\_{i}}{\partial \mathbf{x}\_{j}}\right) = -\mathbf{g}H\frac{\partial \mathcal{L}}{\partial \mathbf{x}\_{i}} + \frac{1}{\rho}\sum\_{j=1}^{2} \frac{\partial}{\partial \mathbf{x}\_{j}} \left(H\tilde{\mathbf{r}}\_{ij}\right) \\ &+ \frac{1}{\rho}\theta\rho\_{a}W\_{i}\left(\left(W\_{1}\right)^{2} + \left(W\_{2}\right)^{2}\right)^{\frac{1}{2}} - \frac{\mathcal{G}}{C^{2}}U\_{i}\left(\left(U\_{1}\right)^{2} + \left(U\_{2}\right)^{2}\right)^{\frac{1}{2}} + \left(-1\right)^{i+1}fHU\_{i+1}, \end{split} \tag{7}$$

where

$$
\tilde{\pi}\_y = \mu\_H \left( \frac{\partial U\_i}{\partial \mathbf{x}\_j} + \frac{\partial U\_j}{\partial \mathbf{x}\_i} \right), \ U\_3 = U\_1. \tag{8}
$$

Note that /2 is often used for the definition of . Here, however, we follow the definition in [11] for convenience.

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):

$$U\_{\wedge} \begin{pmatrix} \mathbf{x}, t \end{pmatrix} = \mathbf{0} \text{ on } \varGamma\_{\mathfrak{e}} \tag{9}$$

#### Tsunami Propagation from the Open Sea to the Coast http://dx.doi.org/10.5772/63814 65

$$\begin{aligned} \mathcal{L}(\mathbf{x}, t) &= \mathcal{L}\_{\Gamma\_{\mathbf{o}}}(\mathbf{x}, t) \\ \sum\_{j=1}^{2} \tilde{\tau}\_{y} n\_{j} &= \mathbf{0} \end{aligned} \qquad \text{on } \Gamma\_{\mathbf{o}}. \tag{10}$$

Here, *ni* is the component in the *xi* direction of the unit outward normal vector on the boundary. On *Γ*o, o (,) is specified. Later, we improve Eq. (10).

Initial conditions:

3

plane is represented by . The horizontal viscosity constant [Ns/m2

0,

( )

t

ò (5)

direction is represented by *Wi*

+

*i*

(8)

+

(9)

direction

, and the

] is

(6)

(7)

<sup>1</sup> , . *<sup>ζ</sup>*

Note that *H* is assumed to be positive in this chapter, whose guarantee is a difficult mathe-

represented by *μ*H, the wind effect coefficient is represented by *θ*, the air density is represent-

Chezy coefficient [m1/2/s] is represented by *C*. Then, the following viscous shallow-water

( ) <sup>2</sup>

*HU*

1 1 2 2 22 1 2 2 1 2 2 1 2 1

> 3 1 , . æ ö ¶ ¶ =+ = ç ÷ ¶ ¶ è ø

Note that /2 is often used for the definition of . Here, however, we follow the definition

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

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

G

boundary conditions on *Γ*c and *Γ*o and the initial conditions at *t* = 0 are as follows [9]:

<sup>c</sup> ( , ) 0 on *U xt <sup>i</sup>* =

*U U x x*

1

r

+ = ¶ ¶ å *<sup>j</sup> j j*

(() ( ) ) (() ( ) ) ( )

*j i U U*

z

<sup>1</sup> 1 ,

*a i i i*

*<sup>g</sup> WW W UU U fHU <sup>C</sup>*

1

= ¶ ¶

*t x*

z

2 2 1 1

æ ö ¶ ¶ ¶¶ ç ÷ + =- + ¶ ¶ ¶¶ è ø

*U U H U gH <sup>H</sup> t x xx*

= =

+ + - + +-

% *<sup>j</sup> <sup>i</sup>*

*ij H*

t m

given initial conditions and boundary conditions.

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

å å % *i i j ij*

*j j j ij*

*i i <sup>h</sup> H h U u dx <sup>H</sup>* =- =

The fixed density of the layer is represented by *ρ* and the stress component in the *x*<sup>i</sup>

z

ed by *ρa*, the wind velocity component speed in the *xi*

matical problem [9].

equations are derived [9]:

qr

r

in [11] for convenience.

where

acting on the *xj*

64 Tsunami

$$U\_{\wedge} \left( \mathbf{x}, \mathbf{0} \right) = U\_{\wedge 0} \left( \mathbf{x} \right), \; \zeta \left( \mathbf{x}, \mathbf{0} \right) = \mathcal{L}\_{0} \left( \mathbf{x} \right). \tag{11}$$

Here, *Ui*0(*x*) and *ζ*0(*x*) are the initial values of *Ui* 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 . For terms without spatial derivatives, the approximation is carried out by using the corresponding step function . The step function is 1 in the barycentric region around the node *k* and is 0 at other places.

$$\left(\frac{\overline{\boldsymbol{\zeta}}^{n+1} - \overline{\boldsymbol{\zeta}}^{n}}{\Delta t}, \overline{\boldsymbol{\phi}}\_{k}\right) + \left(\sum\_{j=1}^{2} \frac{\partial \left(\hat{H}^{n} \hat{U}\_{j}^{n}\right)}{\partial \mathbf{x}\_{j}}, \hat{\phi}\_{k}\right) = \mathbf{0}.\tag{12}$$

(( ) ( ) ) (( ) ( ) ) <sup>1</sup> <sup>2</sup> 1 1 1 <sup>1</sup> <sup>2</sup> 1 1 1 1 1 2 2 2 2 <sup>2</sup> <sup>2</sup> 1 2 2 1 2 <sup>ˆ</sup> ˆ ˆ <sup>ˆ</sup> , , ˆ ˆ <sup>ˆ</sup> <sup>ˆ</sup> <sup>1</sup> ˆ ˆ <sup>ˆ</sup> , , <sup>1</sup> , , *n n n n i i n n i k jk j j n n n n n j i k k H j i jij nn n a i k i U U <sup>U</sup> <sup>H</sup> H U t x U U gH <sup>H</sup> <sup>x</sup> x xx <sup>g</sup> WW W UU U C* f f z f f m r qr f f r + + + = + + + = æ ö æ ö - ¶ æ ö ç ÷ ç ÷ <sup>+</sup> ç ÷ è ø D ¶ è ø è ø æ ö ¶ æ ö æ ö ¶ ¶ ¶ <sup>=</sup> -- + ç ÷ ç ÷ ç ÷ ¶ ç ÷ ¶ ¶¶ è ø è ø è ø æ ö + +- + ç ÷ è ø å å (( ) ) <sup>1</sup> <sup>1</sup> <sup>1</sup> 1 ,. *k i n n i k fH U* f <sup>+</sup> <sup>+</sup> + æ ö ç ÷ è ø + - (13)

In the above, + 1 is determined from Eq. (12). From the obtained value and , , the value , + 1 is determined by Eq. (13). Here, and , are approximate values of (,) and (,) respectively, at the node *k* after *n* time steps, and Δ*t* represents the size of time steps. In addition, the following symbols are also used:

$$\left(\frac{\partial U\_i}{\partial \mathbf{x}\_j}, \mathbf{v}\right) = \int\_{\Omega} \frac{\partial U\_i}{\partial \mathbf{x}\_j} \mathbf{v} \, d\mathbf{x},\tag{14}$$

$$
\hat{U}\_i'' = \sum\_k U\_{i,k}'' \hat{\phi}\_k, \overline{U}\_i'' = \sum\_k U\_{i,k}'' \overline{\phi}\_k. \tag{15}
$$

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.
