OPTIMUM SCHEME OF ADAPTIVE MESH REFINEMENT FOR TSUNAMI SIMULATION

A long wave model with nesting scheme is widely used for tsunami simulation and impact assessment, generally. The conventional nesting scheme, however, has disadvantages for analyzing larger computational mesh size along the coast and needs expensive computation costs. Adaptive Mesh Refinement (AMR) method can consider the interactions of different grids and reduce computational cost together. This study applied the AMR for calculation of long wave propagations, and the di fferent mesh refinement schemes were applied and optimized by a series of numerical experiments. The model was applied to the 2011 Tohoku Earthquake Tsunami, and the model was validated against the post tsunami inundation survey data.


INTRODUCTION
The 2011 Tohoku Earthquake and tsunami gave a catastrophic and tragic earthquake-tsunami disaster for Japan.The rupture area, assumed to be approximately 450 km × 200 km, generated a tsunami that struck Japan over 2000 km along the coast from Hokkaido to Kyushu (e.g.Mori et al., 2011 ;2012).As the spatial scale of the fault caused by the tsunami were several hundred kilometers, the tsunami inundation area ranged more than 500 km.In contrast, the local assessment of tsunami requires several hundred to a few meters to include influence of local bathymetry and structures in coastal area (e.g., to conduct the risk assessment, evacuation planning and real-time forecasting for tsunami warning).The difference of spatial scale to be considered is the order of 10 4 -10 5 .As a result, the tsunami simulation has to cover from the offshore to the nearshore (or inland) area for single computation.A long wave model with nesting scheme is widely used for tsunami simulation, generally.The general nesting scheme divides an entire computation domain into some sub-regions hierarchically and calculates toward high resolution area.The conventional nesting scheme, however, has following some disadvantages for this kind of tsunami simulations.First, the one-way nesting scheme cannot consider influence of local tsunami deformation effects in the fine grid on the coarse grid system.This can sometimes be a problem for the edge wave propagation.Second, the computation cost is expensive because the temporal integration cost is suppressed by the finest grid resolution through the computation.In addition, decision of the appropriate spatial resolution in all domains is required before starting the computation.
In order to resolve these problems with nesting scheme, we apply an Adaptive Mesh Refinement (AMR) method for long wave propagation, which is dynamic grid allocation method for finite difference method.The AMR can consider the interactions of different grids and reduce computational cost before and after main tsunami wave approaches to coasts.The criteria of mesh refinement for tsunamis, however, needs to optimize depending on the governing equations.This study applies the AMR for long wave propagation scheme, and the optimum mesh refinement scheme is examined by the numerical experiments.The application of the model to the 2011 Tohoku Earthquake tsunami is carried out, and the validity of model is checked comparing with survey data.

METHODOLOGY Governing equations
The governing equations of model are linear long wave in order to simplify the code as a first attempt of AMR scheme application.The linear long wave equation for spherical coordinates is used for the governing equations as follows: where η is the surface displacement, P and Q are the vertically integrated horizontal momentum components, (θ, ϕ) are latitude and longitude coordinates, h is water depth, g is the acceleration of gravity, R is the radius of the earth, and f is the Coriolis force.The computation grid is staggered grid, and the equations are discretized with leap-frog method both time and space.The numerical model was corded from scratch for this study.

Mesh refinement algorithm
The spatial resolution is changed dynamically according to some criterion in order to enhance the accuracy of partial differentiation in the AMR method.The AMR method can allocate high resolution grids locally with small number and optimize the numerical solution by balancing the resolution with the numerical accuracy.The AMR method is simply able to implement to the conventional finite difference method, while making the code about mesh refinement is required.
There are the two typical types of the AMR method, one is a cell-base type and another is a block-base type.The block-base type is selected for easiness of making computational code in this study (e.g.Liang. 2010).The spatial resolution of AMR can be determined by the level of mesh at each mesh block.Here ∆ indicates the one side of spatial grid and level means the grid level.As the level increases one, ∆ becomes one half smaller than former grid.
If the starting grid level is one, the relation between the different level in the AMR method have spatial arrangement as shown in Figure 1.
The difference of level between the target mesh and neighbor meshes keeps plus or minus one when the mesh is remeshed to prevent numerical instability of the computation (i.e.2:1 rule).Figure 2 shows an example of mesh locations in different three levels.For example, if the target mesh number is 12 in Figure 2, all neighbor meshes around mesh number 12 should be checked, so that the difference of level between the target mesh and all neighbor meshes is limited plus or minus one.If the difference of level between the neighbor meshes is more than two, the neighbor meshes are refined or coarsed repeatedly until the difference of the level becomes one.The biggest mesh number is limited 2 max to avoid unlimited computational costs.

Spatial interpolation
Conventional finite difference method can be used in the AMR but difference of spatial resolution is necessary to consider for integration.If the level between the neighbor grids is different, the finite difference method is applied based on the interpolated value between the grids.The weighted linear interpolation scheme depends on the locations and levels of neighbor points is applied according to the distance (e.g., Liang, 2010).For example, the upwind quantity for the mesh number 63 which is right side of the mesh number 62 in Figure 2, interpolation can be conducted as follows.
The calculation of the block-base type AMR scheme is simpler than the cell-base type AMR scheme.

Outline of mesh refinement and time-integration
The AMR method is used as a part of the conventional finite different method as shown in Figure 3. Left side shows the computation flow of the conventional temporal integration by finite difference scheme.In the case of consideration the AMR method, we need to additional procedure to the conventional temporal integration as indicated by the right side panel in Figure 3.The AMR method is composed of mesh refinement scheme and interpolation scheme.First of all, the mesh partition is implemented after checking whether or not a target mesh is remeshed.Then, the mesh is recursively optimized by refining the mesh after checking the 2:1 rule.Finally, according to mesh distribution, the interpolation of the wave height, the momentums in x-axis and y-axis and bathymetry are implemented, and the AMR method is inserted before integration scheme.
The mesh refinement is arbitrary depending on the purpose of application and is conducted according to some criterion.Three different criteria for mesh refinement are used in this study.One is surface gradient, the others are ratio of wave amplitude to water depth and Froude number.We examined that whether or not the difference of criteria could affect the computational result through the computation.

RESULTS AND DISCUSSION
The mesh refinement, however, is arbitrary, and optimum scheme for the mesh dividing and merging is necessary for tsunami simulation.A series of numerical experiments of soliton propagation and run-up on a constant slope and hindcast of the 2011 Tohoku Earthquake Tsunami were examined to check the sensitivity of mesh refinement scheme and numerical accuracy.

The computation result with constant slope
In order to check the sensitivity of mesh dividing and merging scheme to numerical accuracy, a series of numerical experiments of soliton propagation on a constant slope was performed.As shown in Figure 4, the input wave profile was one dimensional Gaussian distribution along x-axis as following equation, and where x 0 is the center position of Gaussian distribution, and σ is the standard deviation of the surface.The uniform surface profile was assumed to y-axis although the numerical computation was conducted in two-dimensional domain.The initial velocity/momentum was assumed to be zero for simplicity.Figure 5 shows the maximum water surface elevation over the computational period.The wave was propagated from right side to left side in the figure, and the water surface elevation amplified on constant slope.These computations were executed under some conditions as shown in Table 1.The case 1 is conventional finite difference method without the AMR scheme as a reference test.The mesh level could be increased until two level from the initial condition at most in case2, case3 and case4.That is, if the mesh is divided two steps, the resolution is just the same size as that of case1.As Figure 5 shows, the solution of the wave amplification on the slope approaches theoretical solution, Green's law, from about the point of depth 35 m where the mesh dividing is just started in the three AMR cases.However, we must pay attention to one point as follows.The other mesh dividing and merging criteria were changed the numerical performance, because the mesh dividing or merging is decided only by the threshold.Thus, the accuracy of the AMR can be improved, if the spatial resolution increases by setting the appropriate threshold.For mesh merging, setting the lower threshold less than the upper for mesh dividing can cause decreasing the resolution continuously.On the other hand, total number of mesh is directly gave impact on the computational costs.Therefore, the appropriate mesh dividing and merging criteria depending on the governing equation and the target of application are necessary.This AMR model was applied to tsunami run-up at Kamaishi and Ryoishi bay in Iwate Prefecture by the 2011 Tohoku Earthquake Tsunami to examine the real local tsunami performance practically.The numerical analysis for Kamaishi case was conducted by the conventional finite different method (CFDM) and the AMR method (Figure 6).The computational domain covers two bays and the size is 10 km in the north-south direction and 12 km in the east-west direction, respectively.The bathymetry data is given by 50 m resolution coordinated by the Central Disaster Management Council of Japan.The coastal breakwaters were not resolved due to the limitation of spatial resolution but the offshore breakwater at the bay mouth of Kamaishi is considered correctly.The breakwater at the bay mouth of Kamaishi bay was partially destroyed by the tsunami, but the process of breakdown was excluded for analysis.The detail of offshore breakwater performance at Kamaishi bay was analyzed by Yoneyama et al. (2012).Wave input is the recorded GPS wave gauge data at the offshore of South Iwate by MLIT/PARI (Kawai et al., 2011).The simulated maximum surface elevation was compared with the post-event survey data by the 2011 Tohoku Earthquake Tsunami Joint Survey Group (Mori et al., 2012).The left side number of threshold in Table 2 indicates that the mesh is merged if the A/h of a target mesh is less than this threshold value.The right side number of threshold in the table indicates that the mesh is divided if the value A/h of a target mesh is more than this number.

The 2011 Tohoku Earthquake Tsunami hindcast
Figure 6 illustrates that the survey data marked with circles and the computed maximum wave height by color map.In the case of the AMR computation, we changed the condition of mesh refinement by using the ratio of wave amplitude to water depth M ϵ .The threshold for dividing mesh ranges from A/h=1/5 to A/h=1/20, the threshold for merging meshes ranges from A/h=1/50 to A/h=1/200, respectively.
The major feature of the computed maximum water elevations is accuracy to the observed results.The computational results are basically reasonable in inner part of Kamaishi bay, but it is underestimated around the breakwater.The computational results around the offshore breakwater are not sensitive to the threshold of the mesh refinement.It is considered that the M ϵ in the outside of the bay is not more than the upper threshold for dividing mesh in all cases.The mesh refinement scheme M ϵ is easily able to divide mesh in the coastal area where the depth is small.On the other hand, it is difficult to consider the offshore breakwater effects where the water depth is deep.Thus, the current dividing and merging criteria were not appropriate for this area.In order to divide mesh in the outside of the bay, the upper threshold for dividing mesh must be decreased.However, using the small threshold is not appropriate in the sense of computational cost because the dividing area will be spreaded over the entire domain.
Figure 7 illustrates the snapshot of level distribution in the process of the computation.The corresponding time of Figure 7 is the time the maximum wave amplitude reached around 140 • east longitude (Figure 8).The divided mesh is concentrated on the coastal area for each case, that is, the spatial resolution increases in this area.In this time, there are many divided meshes in the coastal area in all cases, but the maximum wave cannot be caught sufficiently with divided meshes area.If the maximum wave cannot be caught sufficiently due to coarse grid resolution, the wave cannot propagate into inner area of the bay.As a result, it was found that appropriate mesh division was important not only for the coastal area, but also for the outside of the bay where required additional criteria and threshold value of the mesh refinement.
According to the time series of the total number of grid is shown in Figure 9, the number of mesh of case 3 is two times larger than that of case 2. The computing time is directly affected by the number of meshes.It was found that the divided meshes needed to be merged after maximum wave passes through in order to reduce computational cost while keeping the accuracy.The criteria and threshold value of mesh merging is also need to optimize for the AMR method.

CONCLUSIONS
The Adaptive Mesh Refinement code has been developed for tsunami simulation.The mesh refinement criteria and the threshold for tsunami propagation were discussed by two different numerical experiments.The analysis of the 2011 Tohoku Earthquake tsunami was conducted by the conventional finite different method and the AMR method.The AMR method gave lower computational cost than the CFDM.The AMR can show detail of tsunami deformation in nearshore with finer mesh and mesh refinement by amplitude/depth ratio, scheme gives reasonable accuracy.The AMR method could reduce the computational cost keeping accuracy of tsunami wave propagation.However, further study is needed to improve the scheme, upgrade the governing equations into non-linear equations.

Figure 2 :
Figure 2: An example of mesh distribution in the different level grids (the number is the AMR mesh number)

Figure 3 :
Figure 3: Flow chart of computation procedure (left side: main computational scheme, right side: adaptive mesh refinement scheme)

Figure 4 :
Figure 4: Initial condition of run-up on constant slope test