Numerical analysis of impacts of 2011 Japan Tohoku tsunami on China Coast*

REN Zhi-yuan (任智源)

Department of Engineering Mechanics, Shanghai Jiao Tong University, Shanghai 200240, China, E-mail: renzhiyuan815@sina.com

WANG Ben-long (王本龙)

Department of Engineering Mechanics, MOE Key Laboratory of Hydrodynamics, Shanghai Jiao Tong University, Shanghai 200240, China

FAN Ting-ting (范婷婷)

National Marine Environmental Forecasting Center, Beijing 100081, China

LIU Hua (刘桦)

Department of Engineering Mechanics, MOE Key Laboratory of Hydrodynamics, Shanghai Jiao Tong University, Shanghai 200240, China

State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

Numerical analysis of impacts of 2011 Japan Tohoku tsunami on China Coast*

REN Zhi-yuan (任智源)

Department of Engineering Mechanics, Shanghai Jiao Tong University, Shanghai 200240, China, E-mail: renzhiyuan815@sina.com

WANG Ben-long (王本龙)

Department of Engineering Mechanics, MOE Key Laboratory of Hydrodynamics, Shanghai Jiao Tong University, Shanghai 200240, China

FAN Ting-ting (范婷婷)

National Marine Environmental Forecasting Center, Beijing 100081, China

LIU Hua (刘桦)

Department of Engineering Mechanics, MOE Key Laboratory of Hydrodynamics, Shanghai Jiao Tong University, Shanghai 200240, China

State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

(Received October 19, 2012, Revised January 24, 2013)

On the 11th of March, 2011, a subduction earthquake of magnitude Mw9.0 happened at the northeast of Japan, generating a tsunami which resulted in huge damage in Japan. Okada’s elastic fault model is used to generate the deformation of the sea bottom based on USGS sources and UCSB sources respectively. The shallow water equations are solved by the adaptively refined finite volume methods so that it can compute the propagation of tsunami in the Pacific Ocean efficiently. The computed time series of the surface elevation are compared with the measured data from NOAA real-time tsunami monitoring systems for model validation, and UCSB sources derive better results than USGS sources. Furthermore, one nested domain with fine grid and higher topography resolution is combined to compute numerically this tsunami spreading in the Bohai Sea, Yellow Sea, East China Sea, and North of South China Sea. The impacts on China Coast and seas are analyzed and discussed. The results show that the tsunami has almost no impact in the Bohai Sea and Yellow Sea. It has some kind impact on the East China Sea and South China Sea. However, maximum wave height on China Coast is smaller than 0.5 m. It is thus concluded that the 2011 Tohoku tsunami did not generate a significant influence on China Coast.

Tohoku tsunami, numerical simulation, shallow water equation, China Coast

Introduction

On the 11th of March, 2011, an earthquake of magnitude =Mw0.9 happened at the northeast of Japan, generating a tsunami which resulted in huge damage in Japan. For the reason of the high quality of earthquake preparedness in Japan, the earthquake itself caused little damage. However, the subsequent tsunami resulted in the mortality of more than 15 000 people and extensive damaged coastal settlements and infrastructures[1]. Besides, the radiation leak of Fukushima Nuclear Power Plant caused by the tsunami triggered tremendous environment damage which is difficult to estimate[1]. The tsunami also struck the most of the coast of the Pacific rim countries and regions. Therefore, it is vital to simulate the tsunami and predict the impact on China Coast.

Recently, several scholars have simulated this tsunami and analyzed its impacts. Popinet[2]took the Gerris model to simulate the inundation area. This model is based on the shallow water equations. Thenumerical algorithm has some additional capabilities of quadtree-based adaptive discretisation, orthogonal coordinates and parallelism. Zhang et al.[3]used a tsunami model named GeoClaw to simulate the 2011 Tohoku tsunami in the near field and the inundation to the Sendai Airport and Fukushima nuclear power plant with high performance computers with manycore architectures. Grilli et al.[4]developed a new tsunami source by using 3-D finite element model and simulated the tsunami with a Boussinesq model. Saito et al.[5]estimated the initial tsunami water height distribution by the inversion analysis based on the dispersive tsunami simulations. Watanabe et al.[6]calculated the maximum wave heights along Japan Coasts with the shallow water equations.

This article adopts a tsunami model named GeoClaw[7]which was developed by the University of Washington to simulate the 2011 Tohoku tsunami. The model is based on the nonlinear shallow water equations, taking into account the nonlinear effect of the tsunami wave propagation in coastal oceans, whereas it does not involve the dispersion effect. It is based on the tsunami wave height tracking to determine whether to refine mesh. In view of the huge tsunami destruction, rapid tsunami hazard and risk assessments are necessary. Therefore, this tsunami model which ensure the calculation accuracy and calculation efficiency can meet this demand. A numerical investigation of the 2011 Tohoku tsunami based on two earthquake sources (USGS and UCSB) and the assessments of impacts on China coasts are presented, which includes tsunami propagation in the Pacific Ocean for the model validation and tsunami wave height distribution along the coasts of China.

1. Model introduction

1.1 Governing equations

The depth-averaged nonlinear shallow water equations are used to formulate the propagation of tsunamis. In the Cartesian coordinates, the mass conservation equation and the momentum conservation equations can be written as:

in which t is time, h(x,y,t) is total water depth, b(x,y) is the bottom elevation function describing natural bathymetry, (,,)uxyt and (,,)vxyt are the two components of the depth-averaged velocities in the x and y directions. The gravity acceleration is denoted by g. The components of the non-linear bottom friction term are

where n is the Manning coefficient, representing the roughness of the bottom.

Table 1 Parameters of the tsunami fault zone

The shallow water equations can be formulated as the more general form of hyperbolic systems

where q is the vector of unknowns, ()fq and ()gq is the vector of corresponding fluxes, and s is a vector of source terms:

1.2 Numerical aspects

The finite volume method is adopted to discretize Eq.(5) together with the first-order Godunov method

Fig.1 Seafloor displacement (m)

Fig.2 Topography and DART buoys position

Fig.3 The topography and measured points

where Qniis an approximation to the average value of the solution in the i-th grid cell, and both Fni-1/2and Gni-1/2are numerical fluxes approximating the timei -1/2average of the true flux across the left edge of the grid cell over the time interval.

Table 2 Depth of the measured points

The Godunov method uses the Riemann solution to determine cell interface numerical fluxes at each time step. The Riemann solver must also return numerical flux. Thus, we can obtain the form

Table 3 Leading Wave height with and without friction (m)

If the additional terms are introduced, the model can be considered as a second-order algorithm

1.3 Adaptive mesh refinement

The Adaptive Mesh Refinement (AMR) algorithm was developed in the GeoClaw model. This approach was mainly discussed by Berger and Oliger[8], Berger and Colella[9]. The details of implementation of GeoClaw was described in more detail by Leveque et al.[7]

A single coarse (level l) grid comprises the entire domain, while grids at level l+1 are finer than the coarser level l grids with refinement ratios rlxand rly

The essential AMR algorithm can be summarized in the following steps. The model firstly takes a time step of length Δtlon all grids at level l. Using the solution at the beginning and end of this time step, perform space-time interpolation to determine ghost cell values for all level +1l grids at the initial time and all Δtl-1 intermediate times, for any ghost cells that do not lie in adjacent level +1l grids. The ghost cells are defined as additional two rows of grid cells. Take Δtltime steps on all level l+1 grids to bring these grids up to the same advanced time as the level l grids. For the grid cell at level l that is covered by a level +1l grid, replace the solution Q in this cell by an appropriate average of the values from the rlxrlygrid cells on the finer grid that cover this cell. Adjust the coarse cell values adjacent to fine grids to maintain conservation of mass.

2. Fault zone parameters of the tsunami

In the hours following the event, the U. S. Geological Survey (USGS) provided the focal mechanism solutions of the earthquake in terms of its location, size, faulting source, and shaking and slip distributions, in order to place the disaster for timely human response. The Japan Meteorological Agency (JMA) reported the location of the epicenter (38.103oN, 142.861oE). We take the W-Phase moment solution of USGS, which includes strike, dip, rake, depth and the moment M0= 3.9×1022Nm. The length and width of the fault are 400 km and 200 km respectively[10]. Thus, the amount of slip motions can be computed with the equation[11]

where M0is the moment of an earthquake, μ=3.0× 1010N/M, D is the amount of dislocation, L is the length of the fault plane and W is the width of the fault plane. The fault parameters of the tsunami are summarized in Table 1.

Fig.4 A scenario of the tsunami propagation

With the all fault parameters shown in Table 1, the seafloor deformation can be calculated by Okada’s elastic fault model[12]. The results are shown in Fig.1. Figure 1(a) depicts the contours of seafloor displacement, where the solid line denotes uplift, the interval is 1 m, and the dashed line denotes subsidence. Figure 1(c) provides the 3-D view of the seafloor displacement. The maximum displacement is 6.91 m and the minimum is –2.85 m.

Shao et al.[13]provided seismic inversion sources (UCSB) with the JMA hypocenter location, which are based on the analysis of 27 teleseismic broadband P waveforms, 23 broadband SH waveforms, and 53 long period surface waves. The earthquake source plane was divided into 190 25 km by 20 km subfaults (Model phase II). The total seismic moment of UCSB sources is 5.84×1022Nm. Each of 190 subfaults hasbeen calculated with Okada’s model. Figure 1(b) describes the contours of seafloor displacement, which are the linear superposition of the 190 subfaults. And Fig.1(d) provides the 3-D view of the seafloor displacement. The maximum displacement is 15.89 m and the minimum is –4.08 m.

Fig.5 Comparisons of numerical results and DART buoys measurements

3. Simulation of the tsunami and discussion

3.1 Computational domain

The computational domain covers the Pacific Ocean, as shown in Fig.2, whose topography resolution is 2 min. It can be obtained from the National Geophysical Data Center. Additionally, one nested domain for investigating the impacts of the tsunami on China coasts which is from 110oE to 150oE, 20oN to 45oN is considered in the simulation, whose data resolution is about 1 min. In the nested domain, this article takes five levels to simulate the tsunami, whose ratio for AMR is 1:2:4:4:2. The first level grid is set 150× 75, which denotes grid’s size of first level is 1oand the grid’s size of last level is 0.94 min. Four levels are taken in the simulation outside the nested domain for simulating the tsunami propagation in the Pacific Ocean, and the grid size of the last level is 1.88 min. This means the model will take four levels outside the nested domain, and when the tsunami travels into the nested domain, the simulation would take five levels for fine simulation. To control the AMR adaptive refining region, the wave-tolerance that means the difference between the surface elevation and sea level, which decides the AMR regions, is set to be 0.01 m. The simulation will take USGS and UCSB sources respectively.

Real-time surface elevations measured by six DART buoys are obtained from The National Oceanic and Atmospheric Administration (NOAA) real-time tsunami monitoring systems. The topography and dart positions are shown in Fig.2, where the numbers stand for buoy number and the star refers to the epicenter. The border area refers to nested domain which is shown in Fig.3. Ten points near important cities, five locations in main sea regions on China coast and three gauges with numbers (1-Shenjiamen, 2-Shipu, 3-Shantou) are considered in the simulation. The topography and measured points are shown in Fig.3. The depth of the measured points is presented in Table 2.

Fig.6 Contours of arrival time in Pacific Ocean

Fig.7 Comparisons of numerical results and gauge measurements

We carried out numerical simulation of the tsunami propagation with the bottom friction and without the bottom friction respectively. When the friction is considered, the Manning coefficient is 0.014 in the region where water depth is smaller than 1000 m. The comparison of the leading wave height is shown in Table 3. The results show that there is no difference between these two cases. In addition, the wave heights of leading waves based on UCSB are closed to measured data excepted DART_21415, DART_46411, DART_51407 and Shipu. The friction is not involved in the following simulation. It should be noted that the Manning coefficient could have a significant effect on tsunami wave run-up and inundation.

3.2 Propagation across Pacific Ocean

A scenario of the tsunami propagation in the Pacific Ocean with the UCSB sources in 16 h is plotted in Fig.4. The tsunami wave fields are shown at every hour. The tsunami firstly formed in the northeast of Japan, then spread towards the Pacific Ocean. The land boundary has reflection effect in some regions, such as Japan Island, north of Indonesia. The leading wave has diffraction pattern when the tsunami travel across Hawaii.

Comparisons of the computed surface elevation with the measured data at six DART buoys are shown in Fig.5. The dot line refers to data from DART buoys, the dashed line shows the simulation results based on the USGS sources and the solid line denotes numerical results with the UCSB sources. Basically, the results of numerical simulation match reasonably with the real-time data, including the wave height of the leading wave and the arrival time. It reveals that the GeoClaw model can simulate the tsunami amplitude and arrival time effectively. It is interesting to note that the UCSB sources may bring better results than USGS sources in terms of the wave trough, which indicates that exact earthquake source is important for modeling generation of the tsunami. It should be noted that the surface elevation measured at DART_21401 and DART_21415 within the first hour after the earthquake should be caused by the local movement of sea bottom. We did not consider the response of sea bottom outward, and simulate the vertical movement of sea bottom within the epicenter for tsunami generation. Therefore, this kind of oscillation at high frequency cannot be captured in our simulation.

In Fig.6, the contours of arrival time are plotted, where the time interval is 1 h. Numerical results show that, in about 2 h, the tsunami will spread throughout the Japan Island, and arrive at Philippines and Indonesia in 5 h. The tsunami wave will arrive at the west coast of America after about 10 h.

3.3 Impacts on China Coast

When considering the impacts of the tsunami on China coasts, three gauges are selected for comparison. The comparisons of computed results and gauge data are shown in Fig.7. The wave heights of leading wave modeled based on USGS sources are a little larger than the measured data, while the wave heights with UCSB sources match well with measured data.

A scenario of the tsunami’s propagation towards China coasts with the UCSB sources at 16 h after the earthquake is plotted in Fig.8. The tsunami formed in the northeast of Japan, partly travelled towards China. Due to the shadow effects of the Japan Island, the tsunami only could travel to the East China Sea firstly, then to the Yellow Sea and South China Sea. So the impact of the tsunami on the East China Sea is more significant than its effects on other regions. The wave focused at the east of Taiwan Island.

Fig.8 A scenario of the tsunami propagation

Fig.9 Computed time series of wave profiles near the coastal cities

Time series of the surface elevation at the measured points near the important cities, which are based on USGS sources and UCSB sources, are plotted in Fig.9. For the reason that the wave heights near Tianjin and Dalian are around 0.01 m, the time series of the surface elevation are not provided here.

The surface elevations based on the USGS sources are larger. There are minimal impacts on the waters near Tianjin, Dalian and Qingdao, where the maximum wave height is only a few centimeters. The tsunami wave heights near Taiwan do not reach 0.3 m. Due to the barrier effect of the Island of Taiwan, the wave height near the Xiamen sea area does not reach 0.1 m. Maximum wave height near Shanghai is smaller than 0.4 m. The wave heights in Hongkong and Macao are around 0.1 m. It turns out that the impacts of the 2011 Tohoku tsunami on China’s coastal areas is not significant.

Time series of surface elevation in China Sea regions, including the Yellow Sea, East China Sea, Taiwan Strait, and South China Sea are presented in Fig.10. Because the maximum wave height of Bohai Sea is around 0.01 m, it is not included. It is found that the tsunami has little effect on the Bohai Sea and Yellow Sea, whose wave heights are around several centimeters. Because the tsunami needed to travel longer distance due to the terrain. Most energy may dissipate on the way to these seas. The East China Seaendures more impacts, where the wave height reach up to 0.24 m with USGS sources, and the South China Sea and Taiwan Strait’s wave heights are around 0.1 m with the USGS sources.

Fig.10 Computed time series of wave profiles in China Sea regions

The contours of arrival time in 10 h are shown in Fig.11, where the time interval is 1 h. It is interesting to note that the tsunami waves generated in the northeastern Japan spreads and arrives in Taiwan after 4 h, reaches Pearl River Delta areas after 7 h, and then arrives in Shanghai after 10 h.

It appears that the nonlinearity is extremely small and is perhaps negligible. We take one buoy of DART for analysis. The Buoy 21413 is located in the west of the Pacific Ocean and more than 1 000 km from Japan Island with water depth of 5 886 m. The tsunami wave height measured by the buoy is around 0.7 m. The relative water depth is around 1.2×10–4, which refers to small nonlinearity. When tsunami propagates and runs up, the dispersion effect would be important and will be considered in future work.

Fig.11 Contours of arrival time near China Coast

4. Concluding remarks

A numerical investigation of the 2011 Tohoku tsunami’s propagation in the Pacific Ocean and impacts on China with USGS sources and UCSB sources is presented. The numerical model is validated based on the measured time series of the surface elevation from the NOAA DART buoy system, and the UCSB sources can lead to better numerical results than the USGS W-Phase sources. A fine nested grid is used to analyze the impacts of the tsunami on China Coast and sea regions. The results reveal that the tsunami has almost no impact in the Bohai Sea and Yellow Sea. It has some kind impact on the East China Sea and South China Sea, while wave heights along China coasts do not reach 0.5 m. The 2011 Tohoku tsunami did not generate a significant influence on coastal cities in China.

However, both the Okinawa trench in the East China Sea and the Manila Trench in the South China Sea both are potential tsunami sources. It is important to strengthen studies on disaster prevention and mitigation of the local tsunamis for the coastal regions of China, such as the Yangtze River Delta and the Pearl River Delta. Especially, once earthquake occurs in the Manila Trench, the devastating tsunami triggered will attack China coasts and countries around SCS in several hours[14,15], which would result in huge damage.

Acknowledgements

Special thanks go to Prof. Wu Theodore Y. at Caltech, Prof. He You-sheng at Shanghai Jiao Tong University, and Prof. Madsen Per A. at DTU for their encouragement and suggestions on tsunami studies. The authors are grateful to both referees for their con-structive suggestions on revising the paper.

白术中白术内酯Ⅰ、Ⅱ、Ⅲ的靶标预测及蛋白互作网络研究…………………………………………………… 刘志强等（23）：3241

[1] BRUMFIEL G. Fallout forencies hike radiation toll[J]. Nature, 2011, 478(7370): 435-436.

[2] POPINET S. Adaptive modeling of long-distance wave propagation and fine-scale flooding during the Tohoku tsunami[J]. Natural Hazards and Earth System Science, 2012, 12(4): 1213-1227.

[3] ZHANG S., ZHU A. and SONG S. et al. Parallelization of GeoClaw code for modeling geophysical flows with adaptive mesh refinement on many-core systems[C]. 2011 14th IEEE International Conference on Computational Science and Engineering. Dalian, China, 2011.

[4] GRILLI S. T., HARRIS J. C. and BAKHSH T. T. et al. Numerical simulation of the 2011 Tohoku tsunami: Comparison with field observations and sensitivity to model parameters[C]. Proceedings of the 22th International Offshore and Polar Engineering Conference. Rhodes, Greece, 2012.

[5] SAITO T., ITO Y. and INAZU D. et al. Tsunami source of the 2011 Tohoku-Oki earthquake, Japan: Inversion analysis based on dispersive tsunami simulations[J]. Geophysical Research Letters, 2011, 38(7): L00G19.

[6] WATANABE Y., MITOBE Y. and SARUWATARI A. et al. Evolution of the 2011 Tohoku earthquake tsunami on the Pacific coast of Hokkaido[J]. Coastal Engineering Journal, 2012, 54(1): 1250002.

[7] LEVEQUE R. J., GEORGE D. L. and BERGER M. J. Tsunami modeling with adaptively refined finite volume methods[J]. Acta Numerica, 2011, 20(1): 211-289.

[9] BERGER M. J., COLELLA P. Local adaptive mesh refinement for shock hydrodynamics[J]. Journal of Computational Physics, 1989, 82(1): 64-84.

[10] FURUMOTO A. S., TATEHATA H. and MORIOKA C. Japanese tsunami warning system[J]. Science of Tsunami Hazards, 1999, 17(2): 85-105.

[11] LIU P. L.-F., WANG X. and SALISBURY A. J. Tsunami hazard and early warning system in South China Sea[J]. Journal of Asian Earth Sciences, 2009, 36(1): 2-12.

[12] OKADA Y. Surface deformation due to shear and tensile faults in a half-space[J]. Bulletin of the Seismological Society of America, 1985, 75(4): 1135-1154.

[13] SHAO G., LI X. and JI C. et al. Focal mechanism and slip history of the 2011 Mw 9.1 off the Pacific coast of Tohoku earthquake, constrained with teleseismic body and surface waves[J]. Earth Planets Space, 2011, 63(7): 559-564.

[14] WU T.-R., HUANG H.-C. Modeling tsunami hazards from Manila trench to Taiwan[J]. Journal of Asian Earth Science, 2009, 36(1): 21-28.

[15] DAO M. H., TKALICH P. and CHAN E. S. et al. Tsunami propagation scenarios in the South China Sea[J]. Journal of Asian Earth Science, 2009, 36(1): 67-73.

10.1016/S1001-6058(11)60399-6

* Project supported by the National Natural Science Foundation of China (Grant No. 10972138), the Natural Science Foundation of Shanghai Municipality (Grant No. 11ZR1418200) and the Key Doctoral Programme Foundation of Shanghai Municipality (Grant No. B206).

Biography: REN Zhi-yuan (1986-), Male, Ph. D. Candidate

LIU Hua, E-mail: hliu@sjtu.edu.cn