Proceedings of the International Symposium on Engineering Lessons Learned fromthe 2011 Great East Japan Earthquake, March 1-4, 2012, Tokyo, JapanLARGE-SCALE SIMULATION OF
SEISMIC-WAVE PROPAGATION OF
THE 2011 TOHOKU-OKI M9 EARTHQUAKETaro OKAMOTO1, Hiroshi TAKENAKA2, Takeshi NAKAMURA3 and Takayuki AOKI4
1

Assistant Professor, Department of Earth and Planetary Sciences, Tokyo Institute of Technology,
Tokyo, Japan, okamoto.t.ad@m.titech.ac.jp
2
Associate Professor, Department of Earth and Planetary Sciences, Kyushu University,
Fukuoka, Japan, takenaka@geo.kyushu-u.ac.jp
3
Research Scientist, Earthquake and Tsunami Research Project for Disaster Prevention, Japan Agency
for Marine-Earth Science and Technology, Yokohama, Japan, t_nakamura@jamstec.go.jp
4Professor, Global Scientific Information and Computing Center, Tokyo Institute of Technology,
Tokyo, Japan, taoki@gsic.titech.ac.jp
ABSTRACT: We describe a large-scale simulation of the seismic wave propagationfrom the 2011 Tohoku-Oki earthquake. We combine structure models for topography,
sediments, crust and subducting plates to construct a 3D structure model for thenortheastern Japan area. Then we simulate the near-field wave propagation by using a
FDTD method which is accelerated by the use of multi-GPUs (Graphics Processing Unit).
The simulated wave-field exhibits strongly complex pattern reflecting the sourcecomplexity and the heterogeneities around the source region.
Key Words: Great East Japan earthquake, ground motion, finite-difference method,
topography, heterogeneity, GPU computing.

INTRODUCTION
The Tohoku-Oki earthquake on March 11, 2011 (MJMA 9.0; Fig. 1) generated strong shaking reaching
the maximum intensity (seven) on the JMA’s scale and caused devastating tsunamis with run-up
heights exceeding 30 m. 15,845 people were lost, and 3,368 people are still missing (January 27,2012). Such mega-thrust earthquake was not expected to occur along the plate interface off thenortheastern Japan. Thus it is very important to study this event for understanding the geophysicalcondition of the generation of mega-thrust earthquake, the characteristics of the induced strong groundmotions, and the mechanism of the excitation of the large tsunamis.
The ground motion records of this event are quite important data for the quantitative studies on theearthquake source and the induced damages. However, modeling of the ground motions is not a simple
task because of the strong lateral heterogeneity in and around the Japan trench: steeply varying
topography, oceanic water layer, thick sediments, crust with varying thickness and subducting oceanic
plate can all affect the seismic waves radiated from suboceanic earthquakes (e.g., Okamoto 2002,
Nakamura et al. submitted paper). Thus the structural model is an important factor in the study of349

waveform modeling.
The modeling of the ground motion induced by this event is a computational challenge: large
memory size and fast computing devices are required because the huge fault size of the earthquake(about 500 km×200 km) imposes a very large domain size for the simulation. For example, for a
finite-difference domain of 960 km long, 480 km wide and 240 km deep and for a (constant) grid
spacing of 0.15 km, a quite large grid size of 6400×3200×1600 or 33 billion of grid points arenecessary.

Fig. 1 Aftershock distribution of the 2011 Tohoku-Oki earthquake. The epicenters ofearthquakes (MJMA≧4.0) from Mar. 09, 2011 to Apr. 11, 2011 (JST) are shown.
Earthquakes with magnitude larger than or equal to 7 are shown by red stars.
Hypocentral data determined by Japan Meteorological Agency are used.
Therefore, we need to develop numerical methods that can precisely incorporate the effects of theheterogeneous structure including the land-ocean topography. Further, we need to confirm thefeasibility of the methods in the case of large-scale problem: the computation must be done within atolerable time.
Thus, in this paper we first review our 3-D finite-difference time domain (FDTD) method(Okamoto et al. 2010 and 2011) to describe how we implement the land and ocean-bottom topography,
oceanic layer and other lateral heterogeneity. In order to simulate the wave propagation with a large
grid size, we adopt the GPU (graphics processing unit) computing to our finite-difference program.
Next we show the preliminary structure model for the northeastern Japan region, which weconstructed by combining models of sediment, crust and subducting plates. We then present the resultsof the simulation of the wave propagation based on a preliminary source model of the 2011
Tohoku-Oki earthquake.

350

SIMULATION METHOD
3D finite-difference time domain (FDTD) method
We apply the time domain, staggered-grid three-dimensional finite-difference scheme (e.g., Graves
1996). The components of the particle velocity (vi; i=x, y, z) and the stress tensor (ij; i, j=x, y, z) are
the field variables. We use three material parameters (Lamé coefficients and density) as we assume
isotropic and elastic material for the simulation. Thus twelve variables are assigned to a unit cell (i.e.,a heterogeneous formulation: Fig. 2). The finite-difference is fourth-order in space and second-order intime. We apply the absorbing boundary condition (Cerjan et al. 1985) near the side and the bottom
boundaries, and A1 absorbing condition (Clayton and Engquist 1977) at the bottom. A periodic
condition is imposed on the side boundaries.

Fig. 2 A single unit cell of the staggered grid. Open symbols denote the variables thatbelong to the neighboring cells (from Okamoto et al. 2010b).
Unified implementation for land-ocean topography
Recent studies have been revealing the effects of land-ocean topography and oceanic layer on theseismic ground motions due to suboceanic shallow earthquakes: by using the finite-difference
simulations Nakamura et al., (2009) and Nakamura et al. (submitted) showed that the coda waves were
amplified and elongated when the land-ocean topography and oceanic layer were considered in thesimulation. Therefore, such effects need be incorporated in the simulations of seismic wave
propagation and be studied further for better understanding of the ground motions due to shallow
suboceanic earthquakes.
In order to incorporate such structural effects, we have recently proposed an approach to model
structures with both of the land and ocean topography in 3D seismic modeling with thefinite-difference method (Takenaka et al. 2009). The approach unifies the implementation for irregularfree-surface (i.e., the land topography) proposed by Ohminato and Chouet (1996), and that for
irregular liquid-solid interface (e.g., ocean bottom) proposed by Okamoto and Takenaka (2005).
The liquid-solid interface scheme satisfies the boundary conditions between the elastic waves(solid side) and the pressure waves (liquid side) with an accuracy of O(h) where h being the space
increment. For this, (1) we need to put the boundary (i.e. zero rigidity) thorough the grid points of thecomponents of shear stresses and (2) to apply the standard second-order centered equation withaveraged density (i.e., (1/2)(solid + liquid)) to the normal velocity component at the grid points on theboundary. The second condition also imposes second-order operators instead of the fourth-order ones
near (but only near) the liquid-solid boundary. Okamoto and Takenaka (2005) derived these conditions
for the velocity-stress staggered grid finite-difference scheme, and numerically confirmed the validity
by comparing the results with those of a wavenumber integration method (Okamoto and Takenaka

351

1999) in two-dimensional case. Nakamura et al. (2011) extended the liquid-solid implementation to
three-dimensional case and studied the accuracy by numerical experiments.
For the irregular free-surface, the same conditions (1) and (2) must be imposed: in this case theaverage density becomes (1/2)(solid + air) at land and (1/2)(liquid + air) at ocean. In the method of
Takenaka et al. (2009) the material parameters at the center of the cell are first defined. Second, inorder that the above conditions for material parameters be automatically satisfied, the arithmetic
average values of the densities of two adjoining cells are used for effective values on the faces of thecell, and the harmonic average values of the rigidities of four adjoining cells are used for effective
value on the edges. Nakamura et al. (2009) and Nakamura et al. (submitted) applied this unified
implementation. In the present study, we also apply the unified implementation in the FDTD program
for the GPU computing.
GPU computing
The GPU (graphics processing unit) is a remarkable device due to its many-core architecture and highmemory bandwidth (Fig. 3). The GPU delivers extremely high computing performance at a reduced
power and cost compared to conventional central processing units (CPUs). The simulation of seismic
wave propagation is a typical memory-intensive problem which involves a large amount data transfer
between the memory and the arithmetic units, while the number of arithmetic computations isrelatively small. This is the reason we adopt the GPU computing to the simulation of seismic wave
propagation: it can benefit from the high memory-bandwidth of the GPUs. In the Appendix we
describe how we adopt the GPU computing to our FDTD codes. We used NIVIDIA CUDA C(NVIDIA 2010) in programming the FDTD codes.
We use the TSUBAME 2.0 supercomputer in Global Scientific Information and Computing Center,
Tokyo Institute of Technology, for the seismic wave simulation employed in this paper. The
TSUBAME 2.0 is a large-scale GPU cluster equipped with 1408 nodes and 4224 of NVIDIA M2050
GPUs, and has a peak performance of 2.4 PFlops (peta-flops). It is ranked as the world fifth fastest
supercomputer in the November 2011 TOP-500 list (www.top500.org).

Fig. 3 A simplified diagram of a GPU (NVIDIA M2050) used in this study (modified
from Okamoto et al. (2010b)). A single GPU has 14 processors (multiprocessors) and asingle multiprocessor has 32 processor cores. Thus 448 processor cores are integrated ina GPU. The size of the global memory is 3 GB (giga-byte).
3D STRUCTURE MODEL
We compile several detailed structure models in and around the northeastern Japan region to construct
the 3D model for the simulation. Table 1 summarizes the components of the 3D model and references
for them. All these models are defined with several layers divided by irregular, curved interfaces.
Some of the interfaces of the constructed 3D model are shown in Fig. 4. The total number of the layers

352

used in the 3D model is 19 including the water layer as the top layer and the mantle at the base of themode. The model area (red rectangle on the topography map) includes southernmost part of Hokkaido,
almost all of the northeastern Honshu and Kanto area. It also includes some of the Izu islands in thesouthernmost part of the model. The whole Japan trench and outer-rise area are included in the model.
The 3D model constructed here is a preliminary one: the structure models for the land area and forthe oceanic area have been developed separately so that there are differences in the structuralparameters in the same areas between the models. Some interfaces from different models intersect
each other. We have carefully removed such inconsistencies in compiling the models. It is required to
improve the 3D model in the future work by waveform modeling.
Table 1 Components of 3D model and references
Components
Land and ocean topography
Sediments(Japan Seismic Hazard
Information Station)
Land crust
Pacific plate
Philippine Sea plate

Reference
Kisimoto 1999
Fujiwara et al. 2009
Baba et al. 2006
Baba et al. 2006
Nakamura et al. 2010

Fig. 4 Examples of the interfaces assumed in the 3D model. From left to right: theland-ocean topography, the base interface of the tenth layer (deepest sediment), theMoho, and the top interface of the Pacific and the Philippine Sea plates. The red stardenotes the JMA epicenter of the 2011 Tohoku-Oki earthquake. The red rectangle on the
topography map indicates the horizontal area of the FDTD domain used in thesimulation.
SIMULATION OF TOHOKU-OKI EARTHQUAKEIn this section we discuss about an example of large-scale seismic waveform modeling of the 2011
Tohoku-Oki earthquake by using the GPU-accelerated parallel FDTD method.
FDTD parameters
As described in Introduction, a large-scale computing is required to model seismic waves from the2011 Tohoku-Oki earthquake. As summarized in Table 2, for the example shown here we wererequired to use 1000 GPUs (334 nodes) of the TSUBAME 2.0 supercomputer. This is about 24 % of
the entire resources of the TSUBAME 2.0.

353

The maximum frequency in Table 2 is obtained as follows: based on the argument by Moczo et al.
(2000) we use 6×x = 0.9 km as the minimum wavelength. The minimum S-wave velocity is 0.55
km/s for the topmost sediment. Thus the maximum frequency fmax is: fmax = 0.55÷0.9 = 0.61 Hz. The
area of the top sediment is quite limited and for most of the FDTD domain the seismic waves up to
about 1 Hz are correctly modeled.
Table 2 Parameters for the parallel GPU-FDTD simulation
NX, NY (horizontal grid size) and NZ (vertical)
NT (time steps)
x (grid spacing common for all X,Y and Z)
t (time increment)
Number of unit cells
LX, LY, LZ (dimensions of the FDTD domain)
Computed duration (NT×t)
Maximum Frequency
Number of GPUs
Grid size of a single subdomain
Number of subdomains in X, Y, and Z directions
Total memory
Total computation time
Performance (single precision)
Total size of output files

6400×3200×1600
44000
0.15 km0.005 s32.8 billion
960 km×480 km×240 km
220 s
0.61 Hz
1000 (on 334 nodes)
320×320×320
20×10×5
1.4 TB
5768 s33.2 TFlops
1.0 TB

Fig. 5 Slip distribution of the preliminary source model (Okamoto et al. 2011) used forthe simulation of 2011 Tohoku-Oki Earthquake.
Preliminary source model
We use a preliminary source model determined by ourselves (Okamoto et al. 2011) for the example of354

the simulation shown in this paper. The preliminary source model was determined by using a
non-linear waveform inversion method (Okamoto and Takenaka 2009) applied to the teleseismic
broadband body waveforms. We incorporated the effect of the near-source laterally heterogeneous
structure (including the ocean layer and sediments) on the synthetic waveforms by using a
2.5-dimensional finite difference method. The final slip distribution of the preliminary model is shownin Fig. 5. We obtained a heterogeneous rupture process with large slips off Miyagi prefecture near the
JMA epicenter. This feature is similar to that of the slip distribution obtained by Koketsu et al. (2011).
.

Fig. 6 Snapshots of the simulated wave-field. The vertical component of the particle
velocity is visualized with color scale. Red and blue colors denote upward anddownward motions, respectively. From left to right, top to bottom the snapshot at 10 safter the onset of the rupture, at 30 s, 50 s, 70 s, 90 s, and 110 s, respectively, are shown.

355

Fig. 7 Same as Fig. 6 but for snapshots at 130 s, 150 s, 170 s and 190 s.
The source model consists of point sources distributed on the grid points (23×11) on the fault
plane with a spacing of about 20 km. The source terms (i.e., the moment tensors of each point source)
are applied in the 3D FDTD using the stress components (Okamoto 1994; Olsen et al. 1995).
Results of the waveform modelingFig. 6 and 7 show the snapshots of the ground motion (particle velocity) computed for the preliminary
source model of the 2011 Tohoku-Oki earthquake. The ground motions at the free-surface in the landarea and at the solid side of the ocean bottom in the oceanic area are visualized in color scale.
Although not visualized, the propagation of the pressure waves within the ocean water layer is alsomodeled in the simulation (see section “Unified implementation for land-ocean topography”). We areable to observe strong interference due to the complex source model and the heterogeneous structure.
The complexity is especially significant in the later part of the simulation. It is also remarkable that,from about 130 s to 170 s, regions with strong positive motion (red) and negative motion (blue)
emerge near the coast of the Fukushima prefecture. Motions in the regions overwhelm the wave-field
propagated from around the epicenter, and largely distort the pattern of the wave propagation. The
source that provides this feature is the relatively large slip near the Fukushima prefecture retrieved in
the teleseismic waveform inversion result (Fig. 5). Such a feature was also observed in the strong
ground motion records and the strong-motion-generation-areas (SMGAs) were identified below thecoast of Miyagi to Fukushima prefectures (e.g., Kurahashi and Irikura 2011). These characteristics
must be studied further in the future analysis, with improvements in the source model as well as in the
3D structure model.

356

CONCLUSION
In order to simulate the wave propagation from the 2011 Tohoku-Oki earthquake, we constructed apreliminary 3D structure model for the northeastern Japan region by compiling the models for
topography, sediments, crust, and subducting plates. We adopted the GPU computing to accelerate thelarge-scale finite-difference simulation of the seismic wave propagation. By using about 24 % of theresources of the GPU supercomputer, TSUBAME 2.0, at Tokyo Institute of Technology, Japan, we are
now able to simulate the seismic wave propagation from the source of the huge fault plane up to about
1 Hz within a tolerable computation time. The snapshots of simulated wave-field exhibit strongly
complex pattern because of the complex source model, the irregular topography and the shallow
heterogeneities. Such effect must be considered in improving the source model, the 3D structure
model, and in the future quantitative study on the observed ground motion records.
ACKNOWLEDGEMENTS
We are grateful to K. Kisimoto, H. Fujiwara and T. Baba and the Japan Seismic Hazard Information
Station (J-SHIS) in National Research Institute for Earth Science and Disaster Prevention forproviding the structural models. We are also grateful to Global Scientific Information and Computing
Center (GSIC), Tokyo Institute of Technology, for providing opportunity (Grand Challenge) to uselarge number of nodes of TSUBAME. This research was partially supported by Grant-in-Aid for
Scientific Research (KAKENHI 23310122) from Japan Society for the Promotion of Science (JSPS).
APPENDIX: GPU ACCELERATION FOR FDTD SEISMIC WAVE MODELILNG
We here briefly review how we adopt GPU computing to the simulation of seismic wave propagation(see Okamoto et al. 2010 and Okamoto et al. (in press) for detail).
Memory optimization in GPU computing
The GPU is characterized by its hundreds of cores for arithmetic instructions (Fig. 3): we need toprovide sufficient amount of data to these cores. The bandwidth of the main memory (called “global
memory” in NVIDIA GPUs) connected to the GPU is much faster than that of the memory system
connected to the conventional CPUs (e.g., Kirk and Hwu 2010). However, 400 to 600 clock cycles ofmemory latency still occurs in transferring the data between the global memory and the GPU. It is thus
critical for performance to optimize the memory manipulation in memory-intensive applications suchas the simulation of seismic wave propagation. Therefore, we use the fast (but small) memories in theGPU, the registers and the shared memory, as software managed cache memories to reuse the data andhence to reduce the data transfer from the global memory (Abdelkhalek et al. 2009; Micikevicius
2009; Michéa and Komatitsch 2010; Okamoto et al. 2010a,b; see Fig. 8 for the method).
In order to reduce the access to the global memory, we define the material parameters only at thecenter of the unit cell. The material parameters at the grid points for particle velocities and shear
stresses are computed by using the values defined at the center of the unit cells at every time steps.
The unified implementation method mentioned above (Takenaka et al. 2009) allows us to compute theappropriate values of the material parameters. We further use the look-up table method for the material
parameters: instead of the three material parameters (density and Lamé coefficients) only a single
index is required for a unit cell in the look-up table method. This also reduces the access to the global
memory. In addition to these techniques we optimize the number of grid points (i.e., the block size)
assigned to a single execution unit (called a multiprocessor), because the memory transfer rate is betterfor grouped memory transaction using blocks of proper size (typically 16 to 64) than that for serialized,357

one-by-one memory transaction.
Parallel computing with multi-GPUs
For parallel computing with multi-GPUs, we decompose the FDTD domain into subdomains andallocate a single subdomain to a single GPU. This is necessary for large-scale simulations because thesize of the global memory of a GPU is not large (e.g., 3 GB in the case of NVIDIA M2050). We here
use the MPI library for our parallel FDTD program (i.e., we apply a MPI only parallel programming
model).

Fig. 8 Example of the use of the shared memory and registers. Values on red points are
stored in shared memory and those on gray point in registers. In the computations fork-th depth the stress component xz is stored in register. For the next depth the value inthe register is moved to the shared memory to reduce the data transfer from the global
memory (from Okamoto et al. 2010b).

Fig. 9 Schematic illustration of the 3D domain decomposition (from Okamoto et al.
2010b). Blue regions indicate the ghost zones.
In the parallel computing we need to exchange the data in the ghost zones between theneighboring subdomains. We adopt the three-dimensional (3D) domain decomposition (Fig. 9). In thisdecomposition method the whole domain can be extended in all three Cartesian directions, andtherefore the method brings large freedom in defining the shape of the domain. Also, for a fixed wholedomain size, the communication time in the 3D domain decomposition decreases with the increasing
number of the subdomain. This is because the communication time is proportional to the surface area
which decreases with the size of the subdomain. Here we note that, in GPU computing, the data within
the ghost zones in the GPU memory must be once copied to the host memory in order to exchange
them with other GPUs installed on the other nodes. Thus in decomposing the domain, we use
dedicated memory buffers for the ghost zones to improve the data transfer rate between the GPU
memory and the host memory (see Okamoto et al. 2010b and Okamoto et al. in press for detail).
We also overlap the procedures for the communication and the computation to reduce the total
processing time (see e.g., Abdelkhalek et al., 2009; Ogawa et al., 2010). In our program we firstcompute for the grid points at and near the sides of the subdomains for ghost point exchange. Then we358start the exchange procedure for the ghost points and the computations for the internal grid points
simultaneously.
With these optimizations, the performance of our 3D FDTD program scales well with theincreasing number of GPUs: by using 400 nodes of the TSUBAME supercomputer and under thecondition of 2 GPUs activated per a node, a very high single-precision performance of about 50
TFlops is achieved in the case of 800 GPUs (85 % of the complete scalability). Even higher
performance of about 61 TFlops (69 %) is achieved in the case of 1200 GPUs under the condition of 3
GPUs per a node (Okamoto et al. in press). The reduced performance shown in Table 2 (33.2 TFlops)
is primarily due to the overhead of the file output procedure during the computation.
REFERENCES
Abdelkhalek, R., Calandra, H., Coulaud, O., Roman, J. and Latu, G. (2009). “Fast seismic modelingand reverse time migration on a GPU cluster”, International Conference on High Performance
Computing & Simulation, 36-43.
Baba, T., Ito, A., Kaneda, Y., Hayakawa, T. and Furumura, T. (2006). “3-D seismic wave velocity
structures in the Nankai and Japan Trench subduction zones derived from marine seismic surveys”,
Abstr. Japan Geoscience Union Meet., S111-006, Makuhari, Japan.
Cerjan, C., Kosloff, D., Kosloff, R. and Reshef, M. (1985). “A nonreflecting boundary conditions fordiscrete acoustic and elastic wave equations”, Geophysics, Vol. 50, 705-708.
Clayton, R. and Engquist, B. (1977). “Absorbing boundary conditions for acoustic and elastic wave
equations, Bull. Seism. Soc. Am., Vol. 67, 1529-1540.
Fujiwara, H., Kawai, S., Aoi, S., Morikawa, N., Senna, S., Kudo, N., Ooi, M., Hao, K. X.-S.,
Hayakawa, Y., Toyama, N., Matsuyama, H., Iwamoto, K., Suzuki, H. and Liu, Y. (2009). “A study
on subsurface structure model for deep sedimentary layers of Japan for strong-motion evaluation”,
Technical Note of the National Research Institute for Earth Science and Disaster Prevention,
No.337.
Graves, R. W. (1996). “Simulating seismic wave propagation in 3D elastic media using staggered-gridfinite differences”, Bull. Seism. Soc. Am., Vol. 86, 1091-1106.
Kirk, D. B. and Hwu, W. W. (2010). Programming massively parallel processors: a hands-on
approach, 258 pp., Elsevier.
Kisimoto, K. (1999). “Combined bathymetric and topographic mesh data: Japan250m.grd”, Geological
Survey of Japan, Open-file Report, No. 353.
Koketsu, K., Yokota, Y. Nishimura, N., Yagi, Y., Miyazaki, S., Satake, K., Fujii, Y., Miyake, H., Sakai,S., Yamanaka, Y. and Okada, T. (2011). “A unified source model for the 2011 Tohoku earthquake”,
Earth Planet. Sci. Lett., Vol. 310, 480-487.
Komatitsch, D., Michéa, D. and Erlebacher, G. (2009). “Porting a high-order finite-element earthquake
modeling application to NVIDIA graphics cards using CUDA”, J. Parallel Distrib. Comput., Vol.
69, 451-460.
Komatitsch, D., Erlebacher, G., Göddek, D. and Michéa, D. (2010). “High-order finite-element seismic
wave propagation modeling with MPI on a large GPU cluster”, J. Comp. Phys., Vol. 229,7692-714.
Kurahashi, S. and Irikura, K. (2011). “Source model for generating strong ground motions during the2011 off the Pacific coast of Tohoku earthquake”, Earth Planets Space, Vol. 63, 571-576.
Michéa, D. and Komatitsch, D. (2010). “Accelerating a three-dimensional finite-difference wave
propagation code using GPU graphics cards”, Geophys. J. Int., doi: 10.1111/j.1365246X.2010.04616.x.
Micikevicius, P. (2009). “3D finite-difference computation on GPUs using CUDA”, GPGPU-2: Proc.
2nd Workshop on General Purpose Processing on Graphics Processing Units, 79-84, Washington
DC, USA
Moczo, P., Kristek, J. and Halada, L. (2000). “3D fourth-order staggered-grid finite difference359schemes: stability and grid dispersion”, Bull. Seism. Soc. Am., Vol. 90, No. 3, 587-603.
Nakamura, T., Takenaka, H., Okamoto, T. and Kaneda, Y. (2009). “Finite-difference simulation ofstrong motion from a sub-oceanic earthquake: modeling effects of land and ocean-bottom
topographies”, American Geophysical Union, Fall Meeting, S43B-1981, San Francisco, USA.
Nakamura, T., Okamoto, T., Sugioka, H., Ishihara, Y., Ito, A., Obana, K., Kodaira, S., Suetsugu, D.,
Kinoshita, M., Fukao, Y. and Kaneda, Y. (2010). “3D FDM simulation for very-low-frequency
earthquakes off Kii Peninsula”, Abstr. Seism. Soc. Japan, P1-06, Hiroshima, Japan.
Nakamura, T., Takenaka, H., Okamoto, T. and Kaneda, Y. (2011). “A study of the finite difference
solution for 3D seismic wavefields near a fluid-solid interface”, J. Seism. Soc. Japan (Zisin), 2ndSeries, Vol. 63, No. 3, 189-196.
Nakamura, T., Takenaka, H., Okamoto, T. and Kaneda, Y. (submitted). “FDM simulation of seismic
wave propagation for an aftershock of the 2009 Suruga Bay earthquake: Effects of ocean-bottom
topography and seawater layer”, submitted to Bull. Seism. Soc. Am.
NVIDIA Corporation. (2010). “NVIDIA CUDA C Programming Guide”, Version 3.2,
http://www.nvidia.com/.
Ogawa, S., Aoki, T. and Yamanaka, A. (2010). “Multi-GPU scalability of phase-field simulation forphase transition – 5 TFlop/s performance on 40 GPUs”, Trans. IPS Japan, Advanced Computing
Systems, Vol. 3, 67-75.
Ohminato, T. and Chouet, B. A. (1997). “A free-surface boundary condition for including 3D
topography in the finite-difference method”, Bull. Seism. Soc. Am., Vol. 87, No. 2, 494-515.
Okamoto, T. (1994). “Moment tensor in finite-difference calculation”, Programme and Abstracts,
Seism. Soc. Japan, 1994 Fall Meeting, abstract C-06.
Okamoto, T. and Takenaka, T. (1999). “A reflection/transmission matrix formulation forseismoacoustic scattering by an irregular fluid-solid interface”, Geophys. J. Int., Vol. 139, 531-546.
Okamoto, T. (2002). “Full waveform moment tensor inversion by reciprocal finite difference Green’s
function.” Earth Planets Space, Vol. 54, 715-720.
Okamoto, T. and Takenaka, H. (2005). “Fluid-solid boundary implementation in the velocity-stressfinite-difference method”, J. Seism. Soc. Japan (Zisin), 2nd Series, Vol. 57, No. 3, 355-364.
Okamoto, T. and Takenaka, H. (2009). “Waveform inversion for slip distribution of the 2006 Javatsunami earthquake by using 2.5D finite-difference Green’s function”, Earth, Planets and Space,Vol. 61, e17-e20.
Okamoto, T., Takenaka, H. and Nakamura, T. (2010a). “Simulation of seismic wave propagation by
GPU”, Proceedings of symposium on advanced computing systems and infrastructures, 141-142.
Okamoto, T., Takenaka, H., Nakamura, T. and Aoki, T. (2010b). “Accelerating large-scale simulationof seismic wave propagation by multi-GPUs and three-dimensional domain decomposition”, Earth,Planets and Space, Vol. 62, No. 12, 939-942, 2010b.
Okamoto, T., Takenaka, H., Hara, T., Nakamura, T. and Aoki, T. (2011). “Rupture Process And
Waveform Modeling of The 2011 Tohoku-Oki, Magnitude-9 Earthquake”, American GeophysicalUnion, Fall Meeting, U51B-0038, San Francisco, USA.
Okamoto, T., Takenaka, H., Nakamura, T. and Aoki, T. (in press). “GPU-accelerated simulation ofseismic wave propagation”, in GPU Solutions to Multi-scale Problems in Science and Engineering,Yuen, D., Wang, J., Johnsson, L., Chi, C.-H., Shi, Y. (Eds.), 250 pp., Springer.
Olsen, K. B., Archuleta, R. and Matarese, J. R. (1995). “Three-dimensional simulation of a magnitude
7.75 earthquake on the San Andreas fault”, Science, Vol. 270, 1628-1632.
Takenaka, H., Nakamura, T., Okamoto, T. and Kaneda, Y. (2009). “A unified approach implementing
land and ocean-bottom topographies in the staggered-grid finite-difference method for seismic
wave modeling”, Proc. 9th SEGJ Int. Symp., CD-ROM Paper No.37.

360