Tsunami Propagation from the Open Sea to the Coast

Abstract

The tsunami generated by the Great East Japan Earthquake (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 research 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 (1978). First, the viscous shallow-water equation system is derived from the Navier-Stokes equations, based on the assumption of hydrostatic pressure in the direction of gravity. The derived equations have the horizontal viscosity term hereditary from the original Navier-Stokes equations. Next, a numerical finite element scheme is shown. Finally, tsunami simulations of Tohoku-Oki are shown using the above-mentioned approach. Our main concern in this chapter is how to set the boundary condition on the open boundary.

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₁ and x₂ represent directions in the horizontal plane and x₃ represents the vertical direction, and the time [s] is represented by t. Assuming the hydrostatic pressure in the x₃ 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₁ and x₂ directions and the x₃ direction, respectively.

${\displaystyle \sum _{j=1}^3\frac{\partial u_j}{\partial x_j}}=0,$E1

$\frac{\partial u_1}{\partial t}+{\displaystyle \sum _{j=1}^3{u}_j\frac{\partial u_1}{\partial x_j}}=-\frac{1}{\rho }\frac{\partial p}{\partial x_1}+\frac{1}{\rho }{\displaystyle \sum _{j=1}^3\frac{\partial {\tau }_{1j}}{\partial x_j}}+f{u}_2,$E2

$\frac{\partial u_2}{\partial t}+{\displaystyle \sum _{j=1}^3{u}_j\frac{\partial u_2}{\partial x_j}}=-\frac{1}{\rho }\frac{\partial p}{\partial x_2}+\frac{1}{\rho }{\displaystyle \sum _{j=1}^3\frac{\partial {\tau }_{2j}}{\partial x_j}}-f{u}_1,$E3

$0=-\frac{1}{\rho }\frac{\partial p}{\partial x_3}-g.$E4

u(x₁, x₂, x₃, t) represents the fluid velocity component [m/s] in the x_i (i=1−3) direction, p(x₁, x₂, x₃, t) denotes the pressure [N/m²], ρ is the density [kg/m³], τ_ij is the stress component [N/m²] in the x_i direction acting on the x_j plane, f is the Coriolis coefficient [1/s], and g is the acceleration [m/s²] 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).

Figure 1.

Variables of the computational model [9].

Eqs. (1)–(3) [9] are integrated with respect to x₃ 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₁, x₂, t) and the average velocity in the x_i (i = 1, 2) direction is U_i (x₁, x₂, t), then H and U_i are given by the following equations:

$H=\zeta -h,U_i=\frac{1}{H}{\displaystyle {\int }_h^{\zeta }{u}_{i}d{x}_3}.$E5

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 ${\tilde{\tau }}_{ij}$. The horizontal viscosity constant [Ns/m²] 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 x_i direction is represented by W_i, and the Chezy coefficient [m^1/2/s] is represented by C. Then, the following viscous shallow-water equations are derived [9]:

$\frac{\partial \zeta }{\partial t}+{\displaystyle \sum _{j=1}^2\frac{\partial }{\partial x_j}\left(H{U}_j\right)}=0,$E6

$\begin{array}{l}H\left(\frac{\partial U_i}{\partial t}+{\displaystyle \sum _{j=1}^2{U}_j\frac{\partial U_i}{\partial x_j}}\right)=-gH\frac{\partial \zeta }{\partial x_i}+\frac{1}{\rho }{\displaystyle \sum _{j=1}^2\frac{\partial }{\partial x_j}\left(H{\tilde{\tau }}_{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{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}fH{U}_{i+1},\end{array}$E7

where

${\tilde{\tau }}_{ij}={\mu }_H\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right),U_3=U_1.$ E8

Note that ${\tilde{\tau }}_{ij}$ /2 is often used for the definition of ${\tilde{\tau }}_{ij}$. 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₁, x₂) 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_i(x,t)=0\text{on}{\Gamma }_c$E9

$\begin{array}{c}\zeta (x,t)={\zeta }_{{\Gamma }_{\text{o}}}\left(x,t\right)\\ {\displaystyle \sum _{j=1}^2{\tilde{\tau }}_{ij}{n}_j}=0\end{array}\}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{on}{\Gamma }_{\text{o}}.$E10

Here, n_i is the component in the x_i direction of the unit outward normal vector on the boundary. On Γ_o, ${\zeta }_{{\Gamma }_{\text{o}}}(x,t)$ is specified. Later, we improve Eq. (10).

Initial conditions:

$U_i(x,0)=U_{i0}(x),\zeta (x,0)={\zeta }_0(x).$E11

Here, U_i0(x) and ζ₀(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 ${\widehat{\varphi }}_k$. For terms without spatial derivatives, the approximation is carried out by using the corresponding step function ${\overline{\varphi }}_k$. The step function ${\overline{\varphi }}_k$ is 1 in the barycentric region around the node k and is 0 at other places.

$\left(\frac{{\overline{\zeta }}^{n+1}-{\overline{\zeta }}^n}{\Delta t},{\overline{\varphi }}_k\right)+\left({\displaystyle \sum _{j=1}^2\frac{\partial \left({\widehat{H}}^n{\widehat{U}}_j^n\right)}{\partial x_j}},{\widehat{\varphi }}_k\right)=0.$E12

$\begin{array}{l}\left({\overline{H}}^{n+1}\left(\frac{{\overline{U}}_i^{n+1}-{\overline{U}}_i^n}{\Delta t}\right),{\overline{\varphi }}_k\right)+\left({\widehat{H}}^{n+1}{\displaystyle \sum _{j=1}^2{\widehat{U}}_j^n\frac{\partial {\widehat{U}}_i^n}{\partial x_j}},{\widehat{\varphi }}_k\right)\\ =-\left(g{\widehat{H}}^{n+1}\frac{\partial {\widehat{\zeta }}^{n+1}}{\partial x_i},{\widehat{\varphi }}_k\right)-{\displaystyle \sum _{j=1}^2\left(\frac{1}{\rho }{\widehat{H}}^{n+1}{\mu }_H^{}\left(\frac{\partial {\widehat{U}}_i^n}{\partial x_j}+\frac{\partial {\widehat{U}}_j^n}{\partial x_i}\right),\frac{\partial {\widehat{\varphi }}_k}{\partial x_j}\right)}\\ +\left(\frac{1}{\rho }\theta {\rho }_a{W}_i{\left({\left(W_1\right)}^2+{\left(W_2\right)}^2\right)}^{\frac{1}{2}},{\overline{\varphi }}_k\right)-\left(\frac{g}{C^2}{\overline{U}}_i^n{\left({\left({\overline{U}}_1^n\right)}^2+{\left({\overline{U}}_2^n\right)}^2\right)}^{\frac{1}{2}},{\overline{\varphi }}_k\right)+\left({\left(-1\right)}^{i+1}f{\overline{H}}^{n+1}{\overline{U}}_{i+1}^n,{\overline{\varphi }}_k\right).\end{array}$E13

In the above, ${\zeta }_k^{n+1}$ is determined from Eq. (12). From the obtained value and $U_{i,k}^n$, the value $U_{i,k}^{n+1}$ is determined by Eq. (13). Here, ${\zeta }_k^n$ and $U_{i,k}^n$ are approximate values of $\zeta (x,t)$ and $U(x,t)$ 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 x_j},\nu \right)={\displaystyle {\int }_{\Omega }\frac{\partial U_i}{\partial x_j}\nu dx,}$E14

${\widehat{U}}_i^n={\displaystyle \sum _k{U}_{i,k}^n{\widehat{\varphi }}_k},{\overline{U}}_i^n={\displaystyle \sum _k{U}_{i,k}^n{\overline{\varphi }}_k}.$E15

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 ζ₀ (x) = 0 and U_i0 (x) = 0, and the boundary conditions ${\zeta }_{{\Gamma }_{\text{o}}}(x,t)=50$ [m] (0 < t < 60 [s]) at the tsunami-generating area and ${\zeta }_{{\Gamma }_{\text{o}}}(x,t)=0$ [m] at other areas (see 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 2.

A computational model of Tohoku-Oki [9]. (a) Contour map of ζ [9]. (b) Time histories of ζ [9].

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 ζ₀ (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:

$\begin{array}{c}{U}_n(x,t)=c\frac{\sqrt{g\left|h\right|}}{H}\zeta \left(x,t\right)\\ U_t(x,t)=0\end{array}\}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}on{\Gamma }_{\text{o}},$E16

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

The second computational model of Tohoku-Oki.

Figure 4.

Contour map of initial water level.

Figure 5.

(a) Contour map of after 8 min, (b) contour map of after 16 min, (c) contour map of after 24 min, (d) contour map of after 32 min, (e) contour map of after 40 min, and (f) contour map of after 48 min.

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.

Figure 6.

Time histories of ζ.

4. Concluding remarks

Figure 7.

Structural analysis in the case of tsunami attack [9].

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.

Acknowledgments

The first author thanks the JKA RING!RING! Project for the financial support during 2015FY.

Parts of this chapter are reproduced from [9], © 2013 The Author(s); originally published under CC BY license.