2. NUMERICAL SCHEME

2.1 Numerical scheme for linearized equation

For the first step to describe the numerical scheme for the tsunami model, the linearized long wave equation without bottom frictions in one dimensional propagation, Eq.(2.1), is introduced.

(2.1)

Let us introduce the finite difference method to solve the above equation numerically. The finite difference method based upon the Taylor expansion series is shown as follows.

(2.2)

where Δt is the grid interval. We can form the "forward " difference by rearranging Eq.(2.2)

(2.3)

where the first term in the right side of Eq.(2.3) is obviously the finite difference representation for the first order of time derivative at  t=t (see Fig.2.1).

Figure 2.1 Central finite difference representations

The truncation error which have the order of Δt, (O(Δt)) is the difference between the partial derivative and its finite difference representation. Moreover we can rearrange the Taylor expansion series in Eq.(2.2) by replacing Δt  by  +Δt /2  and  -Δt/2  and then we obtain "central" difference with the second order of truncation error.

(2.4)

It is interesting that although the expression of the finite difference representations in Eqs. (2.3) and (2.4) are similar, the order of truncation errors are different. By using the above "central" difference method with the staggered numerical points for water level and discharges, which is called the staggered leap-frog scheme, we can descretize Eq.(2.1) as follows.

(2.5)

For dealing with discrete values in numerical computations, η(x,t) and M(x,t) are expressed for the case of the staggered leap-frog scheme as

(2.6)

where Δx and Δt are the grid sizes in x direction and in time t. The point schematics for the numerical scheme are illustrated in Fig.2.2. The points for water depth h is the same as those for water elevation, η.

Figure 2.2 The point schematics for the numerical scheme

The above finite method provide stable result as long as the C.F.L condition is satisfied:

C (celerity) <  Δx /Δt

Details of the stable condition will be discussed in the Chapter 3.1. Imamura & Goto (1988) investigated the truncation errors in three kinds of typical scheme for long waves simulations and showed that in term of numerical accuracy the staggered leap-frog scheme is the best among them.

2.2 Numerical scheme for convection terms

In the present numerical scheme, an "upwind" difference scheme is applied to the convection terms in order to make the computation stable. The reason why this scheme ensures the stability of computation is explained by taking a simple convection equation in the following:

(2.7)

Here the coefficient C is the propagation velocity and is assumed constant. The arrangement of computation points in the present scheme requires the forward difference scheme for the first order time derivations. This yields

(2.8)

In addition, the central difference is applied to the space derivative.

(2.9)

As a result,  is given by

(2.10)

The solution of Eq.(2.10) is implicitly equivalent the solution of Eq.(2.11) with an truncation error of (Δt2+Δx2).  Substituting Eqs.(2.8) and (2.9) into Eq.(2.7) yields

(2.11)

If the second-order derivative with respect to time is rewritten by using the following relationship (this assumption is valid for the progressive waves),

The solution of Eq.(2.11) is the same as the solution of the following diffusion equation in which the diffusion coefficient is negative.

(2.12)

A negative diffusion works to amplify round-off errors with time leading to  an instability. Therefore, Eq.(2.10) is an unstable difference scheme. The more detail about stable and unstable scheme will be discussed in chapter 3.1.

In order to obtain a stable scheme, the space derivative term is approximated by either forward or backward difference depending on the sign of coefficient C. With the forward difference, we have

and with the backward difference

The corresponding differential equations we are going to solve are within the truncation error of O(Δt2+ Δx2), for the forward difference

(2.13)

and for the backward difference

(2.14)

Therefore, to keep the virtual diffusion coefficient positive (or say to ensure the stability of the computation), we have used the backward difference in case of positive C, and the forward difference in case of negative C, in addition to setting . In other words, the difference should be taken in the direction of the flow.  This is the reason why this scheme is called the "upwind" difference.  Although the leap-frog scheme has the truncation error of the order of Δx2, as long as  the  convection  term  concerns, its order become large as Δx.

2.3 Numerical scheme for bottom friction term

The friction term becomes a source of instability if it is discretized with an explicit scheme. To make the discussion of instability simple, let us consider the following momentum equation without convection terms:

(2.15)

The explicit form of Eq.(2.15) is :

(2.16)

When a velocity become large or a total depth is small in a very shallow water, the absolute of coefficient (amplification factor) of the first term on the right hand side of Eq.(2.16) become more than unity, which leads to numerical instability. In order to overcome this problem, an implicit scheme to set a friction term can be basically introduced. For example, a simple implicit form,

(2.17)

ensures numerical stability, because the amplification factor in Eq.(2.17) is always less than unity. However the effect of friction in shallow water becomes so large that numerical results are dumped. Another implicit form, a combined implicit one to the friction term is given by,

(2.18)

This scheme also gives a stable result. It is, however, noted that the above scheme causes a numerical oscillation at the wave front because the amplification factor could be negative.

We should select the best scheme among some implicit ones to apply the bottom friction term with Manning's roughness. Considering the fact that the numerical scheme of convection terms also involve artificial or numerical dissipation, selection of Eq.(2.17) causes much damping in the result. Therefore the present model uses the combined implicit scheme, Eq.(2.18).