A tsunami simulation of Hakata Bay using the viscous shallow-water equations

The tsunami generated by the Great East Japan Earthquake 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 (Coast Eng Jpn 21:157–171, 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. Next a numerical scheme is shown. Finally, tsunami simulations of Hakata Bay and Tohoku-Oki are shown using the approach.


Introduction
The coastal areas of the 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]. Now, numerical simulations are used to develop disaster-prevention measures to deal with such tsunami disasters. They are used to predict potential future tsunami disasters, to design disaster-prevention facilities such as coastal breakwaters and levees, and in tsunami forecasts to predict tsunami attacks immediately after an earthquake occurs [2,3]. The Central Disaster Prevention Council [4] prepares a basic disaster prevention plan and participates in the determination of important disasterprevention matters. If the Tonankai-Nankai Earthquake occurs, it may cause considerable damage. Consequently, the council has been performing numerical calculations to predict the wave height and the arrival time when the tsunami reaches the coast.
The waveform of a tsunami changes significantly in coastal areas compared with the open sea. This is because the water depth becomes shallow and the wave height becomes high when the tsunami reaches coastal areas. In addition, the wave height is amplified by the complex shapes of the closed-off sections of bays, as well as by the ocean floor topography. The shallow-water long-wave equations are generally used in the numerical simulation of tsunami propagation from the open sea to the coast [5,6]. In this study, we have focused on the viscous shallow-water equations [7] and attempted the computation by the finite element method based on the method proposed by Kanayama and Ohtsuka [8]. Thus far, we have studied the two-layer viscous shallow-water equations [9,10], obtained using the layer model [11], which is a well-known technique for quasi-three-dimensional computations. In these, we showed that the mass conservation law and the energy conservation law are satisfied under physically appropriate conditions [12]. We also found that the two-layer equations correspond to the one-layer equations [7] when the upper and lower layers have the same densities and the upper and lower layer equations are combined. Recently, we have found an interesting survey [13] which deals with sound mathematical topics related to the viscous shallow-water equations.
In this paper, 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 [5]. 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 is given below and simulation results for a Hakata Bay model and a Tohoku-Oki model are presented as computational examples.
In the computational example of Hakata Bay, the tsunami generating area is taken to be the epicenter of the West off Fukuoka Prefecture Earthquake, which occurred in 2005 in the Sea of Genkai located northwest offshore of Fukuoka Prefecture. The validity of this computational method can be confirmed by examining the tsunami propagation velocity based on the numerical analysis results. It can also be confirmed that when the tsunami reaches land in the coastal areas, waves are reflected with islands functioning as breakwaters, suppressing the wave height in our computation.

Derivation of the viscous shallow-water equations
As shown in Fig. 1, we consider the shallow-water long-wave flow in which the wavelength is sufficiently long relative to the water depth [7]. First, the viscous shallowwater equations are derived. Orthogonal coordinates [m] are used, where x 1 and Fig. 1 Variables of the computational model 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 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 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.
Equations (1)-(3) 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 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 i are given by the following equations: The fixed density of the layer is represented by ρ and the stress in the x i direction acting on the x j plane is represented byτ i j . The horizontal viscosity constant [N s/m 2 ] is represented by μ H , the wind effect coefficient is represented by θ , the air density is represented by ρ a , the wind 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 shallowwater equations are derived: Finite element approximation is performed for the viscous shallow-water equations (6) and (7) with given initial conditions and boundary conditions. The computational domain is the two-dimensional 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: Boundary conditions: Here, n i is the component in the x i direction of the unit normal vector on the boundary.
Initial conditions: Here, U i0 (x) and ζ 0 (x) are the initial values of U i and ζ , respectively.

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 . For terms without spatial derivatives, the approximation is carried out by using the corresponding step functionφ k . The step functionφ k is 1 in the barycentric region around the node k and is 0 at other places.
In the above, ζ n+1 k is determined from (12). From the obtained value and U n i,k , the value U n+1 i,k is determined by (13). Here, ζ n k and U n i,k are approximate values of ζ(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: Because of the presence of the basic boundary conditions, (12) holds for all nodes except for the nodes on Γ o and (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 (13) uses Tabata's upwind approximation [8,14]. A mathematical justification for the approximation scheme for linearized equations related to the above scheme was given by Kanayama and Ushijima [15,16] (see also Appendix A for selfcontainedness).

Computational examples
The numerical results for Hakata Bay are shown. The tsunami generating area was taken to be the epicenter of the West off Fukuoka Prefecture Earthquake, which occurred in the Sea of Genkai located northwest offshore of Fukuoka Prefecture on March 20, 2005 (maximum seismic intensity: 6 lower). Because it was a strike slipfault earthquake, no tsunami was observed despite the hypocenter being on the ocean floor. It is therefore noted that our numerical results can not be compared with observed data.
The computation is performed with a Coriolis coefficient f = 2ω sin φ, a plan- . This time step size is found to be too small, because we can get similar results even if a 10 times larger time step size is used for the following mesh data. Details will later be published, referring to the linear stability analysis in Appendix A. Though the full nonlinear stability analysis is very difficult, we believe that the present linear stability condition is useful. In fact, the similar stability condition of [18] in the case of θ = 1 produces a good estimate of the 10 times larger time step size (1 [s]) in the numerical experiments. Regarding the mesh data, the total number of nodes is 11,370, the total number of elements is 22,099, and the ocean floor ordinate is fixed at h = −10 [m]. This constant depth is also not a strong simplification, because the actual depth in Hakata Bay is almost flat. We set the initial conditions ζ 0 (x) = 0 and U i0 (x) = 0, and the boundary conditions  Fig. 2 where the tsunami generating line is shown in different color (red). In general, the tsunami is excited by two ways. The first one is to consider it in the initial condition of the water surface ζ , for which we do not have sufficient input information in this artificial tsunami of Hakata Bay. The second one is to consider it in the boundary condition of ζ as in this paper. In our setting, the computational domain is not so wide that the above approach may be the only way.
Figure 3a-f show the ordinates of the water surface (water level) every 4 [min] from after 4 [min] to after 24 [min]. It is apparent that the tsunami propagates to the closed-off section of the bay. In addition, since inundation of the tsunami is not taken into consideration in our simulation, part of the tsunami is reflected by the island. Figure 4 shows the water level change at the three points A-C in Fig. 5. It is confirmed that the presence of an island near the coast suppresses the wave height on the land side of the island. The distance from the tsunami generation area to the point A is 15 where |h| is the water depth. This is considered to be due to the surrounding land.
Next, the numerical results for Tohoku-Oki are shown. When not specifically defined, the same physical values as for Hakata Bay were used. Regarding the mesh  Since the computational domain is not wide, it looks that there is a reflection from the northern boundary in Fig. 7. This artificial reflection can be removed by suitable boundary conditions on Γ o . Details will later be published. Figure 8 shows the water  Fig. 6. The wave height at the point B is higher than the point A after about 1,500 [s]. When the tsunami reaches coastal areas, the water depth becomes shallow and the wave height becomes high.

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. Using the present Hakata Bay model, in which the tsunami generating area is taken to be the epicenter of the West off Fukuoka Prefecture Earthquake (2005), tsunami propagation from the open sea to the coastal area may be produced. In addition, we have constructed 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 be considered to be a preliminary   Fig. 9. We 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.

A.1 Introduction
This appendix summarizes contents of Kanayama and Ushijima [15,16] for selfcontainedness. First we introduce a linearized system [15] of the viscous shallowwater equations. After the derivation of the model problem, we reformulate the linear initial-boundary value problem within the framework of Hilbert space theory. Then, the existence and uniqueness of the solution of the reformulated problem follows from Hille-Yosida theorem. Next, this appendix concerns the finite element analysis of linearized viscous shallow-water equations [16]. Our scheme is obtained by the use of the "piecewise-liner" continuous finite element approximation for both the waterlevel and the velocity with respect to the space variables, together with the step by step integration method for the time variable. The stability criterion of the scheme and the convergence of the approximate solutions to the exact one are obtained in L 2 sense.

A.2 The derivation of a linearized viscous shallow-water system
Let us consider the following nonlinear stationary problem div(H 0 U 0 ) = 0 in Ω. (17) Supposing that the domain Ω is an open bounded set in the plane R 2 with the C ∞ -class boundary Γ , we set the following boundary condition: We assume the existence of a sufficiently smooth solution pair {U 0 , ζ 0 } of the above problem (16)- (18). Further, let us suppose the existence of two positive constants H 0 and H 0 such that Furthermore we assume U 0 ∈ (C ∞ (Ω)) 2 and H 0 ∈ C ∞ (Ω), where Ω denotes the closure of Ω, and C ∞ (Ω) consists of all infinitely differentiable functions on Ω. In the following, C(Ω) denotes the totality of continuous functions on Ω.
As a model linear initial-boundary value problem, we adopt the following linearized viscous shallow-water equations (20) and (21): with the initial and boundary conditions: whereŨ (x) andζ (x) are sufficiently smooth initial data. In the sequel, conditions (17), (18) and (19) will be assumed, and (16) will not be used.
Introduce the following inner product (u, v) X in (L 2 (Ω)) 3 : Due to (19), (L 2 (Ω)) 3 becomes a real Hilbert space with this inner product ( , ) X , which is demoted by X in the sequel. Let us consider the following X-valued evolution problem: where the linear operator A is defined by with the domain D( A) given by Here, the standard notations of Sobolev spaces are used. Namely, let m be a nonnegative integer. The Sobolev space H m (Ω) of order m on Ω is defined by H m (Ω) = {v : ∂ α v ∈ L 2 (Ω), |α| ≤ m}, where We consider H m (Ω) the Hilbert space with the norm: Let D(Ω) be the totality of infinitely differentiable functions on Ω with compact supports in Ω. The closure of D(Ω) in H m (Ω) is denoted by H m 0 (Ω).

Theorem 1 The linear operator A is a densely defined dissipative operator, that is,
The range R( A−λ) is dense in X. Hence, the closure A of A exists and is a maximal dissipative operator in X, since A is dissipative by Theorem 1. A generates a contraction semi-group of linear operators {e t A } t≥0 acting in X. Therefore, we arrive at the following theorem [15].

Theorem 2
Consider the X-valued evolution problem : a ∈ D( A). (E) Then, (E) has a unique solution u(t) = e t A a satisfying u(t) X ≤ a X (t ≥ 0).

Remark 1
For any a ∈ X, the function u(t) = e tA a is called the generalized solution of the problem (E).

A.4 A finite element semi-discrete approximation
For simplicity of notations, we introduce bilinear forms a(U, V ), b(U, V ), c(ζ, V ) and d(U, η) for U, V ∈ (H 1 0 (Ω)) 2 and ζ, η ∈ H 1 (Ω) as follows: Define the function space V and the bilinear form a(u, v) for u, v ∈ V as follows: Set the following weak formulation (π ) of (E): a(u(t), v) for any v ∈ V and t > 0, u(0) = a ∈ V .
(π ) Let h be a "triangulation" of the domain Ω with the representative length h. The curved elements due to Zlamal [17] are adopted near the boundary Γ , while the triangular elements are used in interior of Ω. We assume that the family of "triangulation" { h : 0 < h ≤ h < ∞} is regular, and satisfies the inverse assumption. Namely, there are positive constants Σ and Λ independent of h such that where h T and ρ T denote the diameters of an element T and the inscribed circle of T , respectively. Let W h be the "piecewise-linear" continuous finite element space constructed by More precisely, η h ∈ W h is a polynomial with degree at most 1 on each triangular element of h , while on a curved element T, η h (F T (x)) is a polynomial with degree at most 1 on a reference triangleT , where F T :T x → x = F T (x) ∈ T is the isomorphism fromT onto T determined by Zlamal's method [17]. Now, define the vector-valued finite element space V h by The following problem (π h ) is said to be the semi-discrete Galerkin approximation of (π ): The set V h , being considered as a closed subspace of X, is denoted by X h , which is a Hilbert space with the inner product ( , ) X . Let A h be the bounded linear operator acting on X h , defined by the relation: Since we have (3.5)-(3.9) in [15] for u = U ζ ∈ V , it holds that This estimate together with the finite dimensionality of X h implies that the operator A h is maximal dissipative. Hence, the exponential function {e t A h : t ≥ 0} forms a contraction semi-group acting on X h . The problem (π h ) is equivalently represented in the following X h -valued evolution equation (E h ): The unique solution of (E h ) is given by u h (t) = e t A h a h , which satisfies The main conclusion of this section [15] is:

Theorem 3 Let u(t)
and u h (t) be the generalized solution of (E) and the solution of (E h ), respectively. If lim h→0 a h − a X = 0, then lim h→0 u h (t) − u(t) X = 0 uniformly in t ∈ [0, T] for any finite T.

A.5 A class of finite difference schemes in time
First, several notations of spaces and operators are prepared to state our schemes. Let X be the Hilbert space (L 2 (Ω)) 2 with the inner product: and let Z be the space L 2 (Ω) with the inner product: The spaces V h and W h are denoted by X h and Z h when they are considered as the subspaces of X and Z , respectively. The product Hilbert space Now our finite difference schemes are given in the following : where τ is a given positive time increment, θ is a parameter satisfying and the initial value u 0 = U 0 ζ 0 is given in X h . It is noted that Ushijima [18] deals with the case of θ = 1, while Kanayama-Ushijima [19] slightly modifies the case of θ = 1 2 . Introduce a i (U, V ) and b i (U, V ) (i = 0, 1) as follows: where ν is equal to μ H /ρ in (8). It is noted thatτ i j in (8) is simply written as τ in Appendix A. Then we have It is easy to see that there are positive constants α, μ, β 0 , β 1 and γ with which the following estimates from (38) to (42) hold true: b 0 (U, V ) ≤ β 0 a 0 (U, U ) 1/2 V X for U ∈ (H 1 0 (Ω)) 2 , V ∈ (L 2 (Ω)) 2 , (40) V X ≤ β 1 a 0 (V, V ) 1/2 for V ∈ (H 1 0 (Ω)) 2 .
Denote u h (0) by a h . Let u(t) be the generalized solution of (E). Under these settings, we have that and that if lim h→0 a h − a X = 0, uniformly in t ∈ [0, T] for any finite T.
Remark 2 Finally we remark that the function u h (t) in Theorem 4 satisfies the error estimate: u h (t) − u(t) X ≤ Ch for any t ∈ [0, T] with an h-independent constant C provided that the solution u(t) is sufficiently smooth, and that the initial value a h is chosen so as to satisfy a h − a X = O(h).
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.