Accurate and efficient modeling of global seismic wave propagation for an attenuative Earth model including the center

Physics of the Earth and Planetary Interiors 200–201 (2012) 45–55

Contents lists available at SciVerse ScienceDirect

Physics of the Earth and Planetary Interiors journal homepage: www.elsevier.com/locate/pepi

Accurate and efficient modeling of global seismic wave propagation for an attenuative Earth model including the center Genti Toyokuni a,⇑, Hiroshi Takenaka b a b

National Institute of Polar Research, 10-3, Midoricho, Tachikawa, Tokyo 190-8518, Japan Department of Earth and Planetary Sciences, Faculty of Sciences, Kyushu University, 6-10-1 Hakozaki, Higashi-ku, Fukuoka 812-8581, Japan

a r t i c l e

i n f o

Article history: Received 6 May 2011 Received in revised form 8 March 2012 Accepted 16 March 2012 Available online 28 March 2012 Edited by George Helffrich Keywords: Global seismology Waveform modeling Axisymmetric modeling Finite-difference method (FDM) Anelastic attenuation Memory variables Singularity Earth’s center

a b s t r a c t We propose a method for modeling global seismic wave propagation through an attenuative Earth model including the center. This method enables accurate and efficient computations since it is based on the 2.5-D approach, which solves wave equations only on a 2-D cross section of the whole Earth and can correctly model 3-D geometrical spreading. We extend a numerical scheme for the elastic waves in spherical coordinates using the finite-difference method (FDM), to solve the viscoelastodynamic equation. For computation of realistic seismic wave propagation, incorporation of anelastic attenuation is crucial. Since the nature of Earth material is both elastic solid and viscous fluid, we should solve stress–strain relations of viscoelastic material, including attenuative structures. These relations represent the stress as a convolution integral in time, which has had difficulty treating viscoelasticity in time-domain computation such as the FDM. However, we now have a method using so-called memory variables, invented in the 1980s, followed by improvements in Cartesian coordinates. Arbitrary values of the quality factor (Q) can be incorporated into the wave equation via an array of Zener bodies. We also introduce the multi-domain, an FD grid of several layers with different grid spacings, into our FDM scheme. This allows wider lateral grid spacings with depth, so as not to perturb the FD stability criterion around the Earth center. In addition, we propose a technique to avoid the singularity problem of the wave equation in spherical coordinates at the Earth center. We develop a scheme to calculate wavefield variables on this point, based on linear interpolation for the velocity–stress, staggered-grid FDM. This scheme is validated through a comparison of synthetic seismograms with those obtained by the Direct Solution Method for a spherically symmetric Earth model, showing excellent accuracy for our FDM scheme. As a numerical example, we apply the method to simulate seismic waves affected by hemispherical variations of P-wavespeed and attenuation in the top 300 km of the inner core. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Coupled to advances in computer architecture, numerical methods have progressed swiftly, and continue to do so every year. In particular, progress in simulation techniques for seismic wave propagation has greatly improved our knowledge of seismic wavefields and Earth’s inner structure. One example is inverting seismic waveform data. This is usually used to determine Earth structure at a variety of spatial scales and requires iterative, linearized inversion with the structural model determined at each iteration, through calculation of synthetic seismograms (e.g., Geller and Hara, 1993). Several numerical methods for seismic waveform modeling have been proposed. However, domain methods such as the finite-difference method (FDM), finite-element method (FEM), ⇑ Corresponding author. E-mail address: [email protected] (G. Toyokuni). 0031-9201/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.pepi.2012.03.010

pseudospectral method (PSM), and spectral-element method (SEM) offer the best ways to model realistic seismograms observed after propagation through complicated heterogeneous media. The domain methods discretize a computational domain into grid points with individual material parameters, then solve discretized wave equations to calculate wave motion at each grid point (Takenaka et al., 1998). We can classify the domain methods by focusing on the dimension of medium assigned to, and the wave dimension calculated by, the computation. Three-dimensional (3-D) modeling calculates 3-D seismic wavefields propagating through a 3-D medium, while two-dimensional (2-D) modeling deals with 2-D wavefields on a 2-D slice of the medium. Recently, an alternative has been used frequently, 2.5-D modeling. This assumes invariance of structure in a direction, and computes 3-D seismic wavefields on a 2-D structural slice including source and receiver. The 3-D modeling produces the most realistic results, taking full account of the distribution of structural heterogeneity, pulse shape,

46

G. Toyokuni, H. Takenaka / Physics of the Earth and Planetary Interiors 200–201 (2012) 45–55

and geometrical spreading effects in 3-D. For example, the SEM is now known as the most prevalent method for computations of full complex 3-D global wave propagation (e.g., Komatitsch and Tromp, 2002a,b; Chaljub et al., 2003). Nevertheless, because of its severe demand on computational resources (computation time and memory), full 3-D modeling of the global seismic wavefield has many restrictions. It is therefore not realistic without a large supercomputing system, such as the Earth Simulator (ES) in Japan. For example, Tsuboi et al. (2008) calculated realistic global seismograms with the spectral-element method that are accurate to 3.5 s and longer, using 4056 processors (507 out of 640 nodes) of the ES. Such computation has great potential in various seismological fields, but is still well beyond normal PC and workstation capacities, even for one computation. In investigations of Earth’s inner structure via waveform inversion, iterative forward computations of the seismic wavefield are required, which forces the use of more convenient methods based on 2-D and 2.5-D modeling. The 2-D modeling provides intuitive realization that ignores the out-ofplane direction, since the majority of seismic energy travels in the plane with source and receiver. In spite of its reduced computation time and memory through consideration of structural parameters and wave motion only on the cross-section, 2-D modeling has several serious disadvantages. For example, it cannot correctly model a point source, 3-D pulse shape, or 3-D geometrical spreading, in addition to its out-of-plane motion underestimate. These drawbacks prevent direct comparison of waveforms between observed and synthetic seismograms produced by 2-D modeling. For the global modeling, in contrast to the drawbacks of the 3-D and 2-D approaches, 2.5-D modeling can correctly treat the 3-D seismic wavefield by accounting for out-of-plane motion, with a small computation time and memory comparable to 2-D modeling by focusing only on a cross-section of the 3-D sphere with rotational symmetry (Toyokuni et al., 2005; Toyokuni et al., 2012). Since it assumes global structures to be axisymmetric with respect to the axis through a seismic source, 2.5-D modeling in spherical coordinates has traditionally been called axisymmetric modeling. As discussed later, traditional axisymmetric modeling has been expanded to treat arbitrary asymmetric structural models within its framework. Consequently, we henceforth call this method spherical 2.5-D modeling, to avoid confusion with the term ‘‘axisymmetric’’. Because of its accuracy and efficiency, spherical 2.5-D modeling based on the FDM had been developed and used in many studies of seismic wave propagation (Alterman et al., 1970; Igel and Weber, 1995, 1996; Chaljub and Tarantola, 1997; Igel and Gudmundsson, 1997; Thomas et al., 2000). Such a moderate sized computation method plays an important role for global seismic waveform modeling, and continuous improvement has still been placed on it. For example, Jahnke et al. (2008) extended the scheme of Igel and Weber (1995) for use on parallel computers with distributed memory architecture. They calculated synthetic seismograms at dominant periods down to 2.5 s for global mantle models, using high performance computers and PC networks. This scheme was used by Thorne et al. (2007) to model SH-wave propagation through cross-sections of laterally varying, lower mantle models under the Cocos Plate. Similar approach can also be seen in works on the SEM (e.g., Nissen-Meyer et al., 2007, 2008). From the practical viewpoint, the FDM can provide more user-friendly schemes since coding and gridding are easy to deal with. Toyokuni et al. (2005) extended conventional spherical 2.5-D modeling to an arbitrary heterogeneous structural model, adopting the newly defined ‘‘quasi-spherical domain’’ instead of the usual spherical domain. Toyokuni and Takenaka (2006a) developed a scheme to implement an arbitrary moment tensor point source in the spherical 2.5-D FDM, using the Fourier expansion of

wavefield variables in the transverse direction. As a numerical example, they simulated wave propagation for an asymmetric structure including a stagnant slab, radiated by a point source with the mechanism of the 1994 deep Bolivia earthquake. Toyokuni and Takenaka (2006b) also used this scheme to simulate the seismic wavefield from the 1994 deep Fiji earthquake, to predict waveforms observed at Antarctica that were propagated through an anomaly of the density and seismic wavespeeds below New Zealand. However, all previous works based on the spherical 2.5-D FDM treated only the elastodynamic equation. In addition, the Earth center had been incorporated only by the acoustic wave simulation (Thomas et al., 2000) because of difficulty on solving singularity at this point. In this paper, we extend the scheme described in Toyokuni et al. (2005) and Toyokuni and Takenaka (2006a), to treat attenuative structures and the Earth center. Our scheme has the following superiority: (1) It enables further realistic computations of global synthetic seismograms retaining small computation requirement of traditional axisymmetric approach; (2) It provides 3-D seismic wavefields which can be directly compared with observations; (3) It uses the FDM which can deal with steep structural change, and provides user-friendly program; (4) It simultaneously calculates complete wavefields at all spatial grid points from the Earth center to the free surface, so we can intuitively understand the wave propagation from a wide view. Because of its small memory requirement, this scheme is also suitable for multi-GPU computations which accelerates large scale FD simulation (e.g., Komatitsch et al., 2009; Okamoto et al., 2010). 2. Methods This section describes our FDM scheme developed for an attenuative Earth model including its center. We introduce anelastic attenuation into the 2.5-D FDM scheme in spherical coordinates ðr; h; /Þ. The anelastic attenuation of seismic waves originating from the deviation of Earth material from a perfect elastic solid can be considered using the viscoelastodynamic equation. Computation of the constitutive viscoelastic relation in the time domain was made difficult because the relation becomes a convolution integral, which requires all history of strain prior to the time in question. Development of schemes with so-called memory variables and their continuous improvement made the FD computation practicable in Cartesian coordinates. This section applies this knowledge to computations in spherical coordinates. Furthermore, we face two problems in the treatment of seismic wave propagation through the center of the sphere: reduced lateral grid spacings around the center that perturb the stability criterion of the FD computation, and the singularity of wave equations in spherical coordinates exactly on this point, because r ¼ 0. For the first problem, we introduce a multi-domain technique, in which several domains consisting of FD grids with different lateral grid spacings are connected in the r direction, forming coarser lateral grids around the center. For the second problem, we develop a scheme using linear interpolation of wavefield variables in the r direction, giving values of particle velocity and stress at the center. 2.1. Incorporation of anelastic attenuation Earth material has both elastic and anelastic properties. The anelasticity of the material is usually represented by the dimensionless Q factor. Materials with a higher Q factor reduce the loss of seismic energy as a wave propagates through them. Since the anelastic property significantly changes the amplitudes and phases of observed seismic waveforms, computation of waveforms should

47

G. Toyokuni, H. Takenaka / Physics of the Earth and Planetary Interiors 200–201 (2012) 45–55

be done with consideration of anelasticity for direct comparison of synthetic seismograms with observations. Such a property of Earth material is mathematically described as a combination of elastic solid and viscous fluid (Liu et al., 1976), so we need to solve the constitutive relation of viscoelastic medium for incorporation of the anelastic attenuation. Such relations in the time domain contain a convolution integral of viscoelastic modulus and strain, which states that the computation of stress at a particular time requires all history of strain before that time. This imposes a huge memory requirement and makes it very difficult to implement anelastic attenuation into a time-domain computation such as the FDM. In the 1980s, however, modulus approximation, which replaces the frequency-dependent viscoelastic modulus by a n-th order rational approximation, offered a breakthrough in the treatment of this phenomenon in the time domain (see Supplementary information). This was brought about by use of appropriate rheological models. These models are used to approximate the Q factor in the frequency domain, within a specified frequency band. Through this process, we obtain n ordinary differential equations for n additional internal variables, which replace the convolution integral. Emmerich and Korn (1987) elaborated the method, adopting the generalized Maxwell body (GMB) as a rheological model, which consists of a finite number of Maxwell bodies and one spring, connected in parallel. Carcione et al. (1988a,b) independently developed a theory for the generalized Zener body (GZB), i.e., a parallel connection of a finite number of Zener bodies (standard linear solids; see Supplementary Fig. S1). Although the GMB and GZB had long been used separately, forming two major styles of simulation with anelastic attenuation, Moczo and Kristek (2005) showed that the rheology of these two models is identical. Reviewing the history of global waveform modeling with anelastic attenuation, we should mention a monumental work by Liu et al. (1976), which estimated implications of velocity dispersion attributable to anelasticity for seismic data and mantle composition, introducing the GZB. Since the work of Carcione et al. (1988a,b) is an application of Liu et al. (1976) to the time-domain FD computation, and such strategy is also used to analyze the attenuation structure of Earth’s inner core (Li and Cormier, 2002) that will be addressed in Section 3.2, we henceforth treat only the GZB. For the GZB, a relation between the Q factor and relaxation times becomes

XL l¼1

Q ðxÞ ¼ X

L l¼1

1 þ x2 sl slr 1 þ x2 ðslr Þ2 ; xðsl  slr Þ

Before moving on to a description of the viscoelastic equation in spherical coordinates, we write the governing equation in Cartesian coordinates following JafarGandomi and Takenaka (2007), to show an outline of our computation in the time domain. In isotropic viscoelastic material, seismic wave propagation in Cartesian coordinates is governed by the following equations.

q

@ v i @ rij þ fi ; ¼ @t @xj

  L @ rij @v k @v j @v i 1P _ ij ; þ dij þ M þ rl  M ¼ ðP  2MÞ L l¼1 ij @t @xk @xi @xj    @r lij 1 @v k @v j @v i ; dij þ Dll þ ¼  l r lij þ ðDpl  2Dll Þ @t @xk @xi @xj sr

1 þ x2 ðslr Þ2

where x is the frequency, L is the number of Zener bodies, and sl and slr are strain and stress relaxation times of the l-th Zener body, respectively (see Supplementary information). These relaxation times carry anelasticity on time-domain simulation of the viscoelastic seismic wavefield, instead of using Q. Since the Q factor is defined for both P- and S-waves, we should obtain appropriate values of the stress and strain relaxation times for both waves, through optimization of Eq. (1) with respect to our desired frequency range. Although the optimization can treat any frequency-dependent Q factors, there is a simple optimization approach for frequency-independent (constant) Q, derived by Blanch et al. (1995), called the s-method. It is well known that Q in the Earth is nearly constant over the seismic frequency range, since Earth material is composed of different minerals, and the attenuation in each is contributed by several physical processes (Moczo et al., 2004). Numerical examples in this paper are also only performed for constant Q (see Supplementary information). For simplification, the same stress relaxation Sl times slr can be used for both wave types slr ¼ sPl r ¼ sr (Robertsson et al., 1994), whereas the strain relaxation times are determined for Sl each wave, yielding sPl  and s .

ð3Þ ð4Þ

where i; j; k 2 f1; 2; 3g; t is time, dij is the Kronecker delta, q is the mass density, P and M are the unrelaxed moduli, Dpl and Dll are the differences between unrelaxed and relaxed moduli for each single relaxation, r lij are the memory variables for the l-th Zener body, and v i ; fi ; rij , _ ij are components of the particle velocity vector, body force and M vector, stress tensor, and first-order time derivatives of the components of the source moment tensor, respectively. In Eqs. (2)–(4), the summation convention has been applied. When we set the relaxed moduli as pR and lR , the unrelaxed ones, P and M, are represented as

P ¼ pR þ dP; M ¼ lR þ dM;

ð5Þ

where dP and dM are the relaxation of moduli that have forms

dP ¼

L L 1P 1P Dpl ; dM ¼ Dll ; L l¼1 L l¼1

ð6Þ

with

Dpl ¼ pR









sPl sSl  1 ; Dll ¼ lR l  1 : l sr sr

ð7Þ

As we described above, the superscripts P and S on the strain relaxation time in Eq. (7) indicate those for P and S waves, respectively. For determination of relaxed moduli, we should use P- and S-wave phase velocities V P0 and V S0 at a reference angular frequency x0 , as

2

pR ð1Þ

ð2Þ

3

7 6 L 7 6 ¼ qV 2P0 R6  7; Pl 4PL 1 þ ix0 s 5 l¼1 1 þ ix0 slr 3 2 7 L 7   7; 4PL 1 þ ix0 sSl 5 l¼1 1 þ ix0 slr 6

lR ¼ qV 2S0 R6 6

ð8Þ

ð9Þ

where i is the imaginary unit. 2.2. Viscoelastic wave equation in spherical coordinates Here, we address the viscoelastodynamic equation for spherical 2.5-D modeling. We extend the 2.5-D elastodynamic equation in spherical coordinates for an arbitrary moment-tensor point source (Toyokuni and Takenaka, 2006a) to treat the viscoelasticity. The /dependence of the moment-tensor sources can be considered through the Fourier expansion of wavefield variables in this direction; the seismic wavefield radiated from an arbitrary-oriented dislocation source is decomposed into five different contributions from five moment-tensor elements, i.e., one axisymmetric source, two vertical dip-slip sources, and two strike-slip sources. Since the equations of motion have the same forms as in the elastic case, as evident in Eq. (2), we describe relations among the stress, particle velocity

48

G. Toyokuni, H. Takenaka / Physics of the Earth and Planetary Interiors 200–201 (2012) 45–55

and the memory variable, which correspond to Eq. (3) in Cartesian coordinates, as follows.

^m @r @ v^ m P  2M @ v^ m ðP  2MÞfm m 2ðP  2MÞ m rrA hA v^ /B þ v^ rA ¼ P rA þ þ r r sin h r @t @r @h L P  2M m 1P ^_ m þ v^ hA cot h þ ^rlm M ð10Þ rrA ; r L l¼1 rrA ^m @r @ v^ m P @ v^ m ðP  2MÞfm m 2ðP  MÞ m hhA hA v^ /B þ v^ rA ¼ ðP  2MÞ rA þ þ r sin h r @t @r r @h L ðP  2MÞ m 1P ^_ m þ v^ hA cot h þ ^rlm M ð11Þ hhA ; r L l¼1 hhA ^m @r @ v^ m P  2M @ v^ m Pfm m 2ðP  MÞ m //A hA v^ þ v^ rA ¼ ðP  2MÞ rA þ þ r r sin h /B r @t @r @h L P 1P ^_ m ^r lm  M þ v^ m ð12Þ hA cot h þ //A ; L l¼1 //A r   L ^m @r @ v^ m 1 @ v^ m 1 1P ^_ m rhA hA rA ^r lm  M ð13Þ ¼M þ  v^ m hA þ rhA ; r @h r L l¼1 rhA @t @r   L ^m @r M @ v^ m fm m 1P ^_ m h/B /B v^ hA  v^ m/B cot h þ ^rlm M ð14Þ ¼  h/B ; r sin h L l¼1 h/B @t @h   L ^m @r @ v^ m fm m 1 m 1P ^_ m r/B /B v^ rA  v^ /B þ ^rlm M ð15Þ ¼M  r/B ; r sin h r L l¼1 r/B @t @r where A; B 2 fC; Sg, with f ¼ 1 when A ¼ C; B ¼ S and f ¼ 1 when A ¼ S; B ¼ C; m 2 f0; 1; 2g is the Fourier expansion order, and ^ m ^0 ^m ^m ^0 ^m ^m ^l0 ^lm ^lm the sets fv^ 0i ; v^ m iC ; v iS g; ff i ; f iC ; f iS g; frij ; rijC ; rijS g; fr ij ; r ijC ; r ijS g, and ^_ 0 ^_ m ^_ m fM ij ; M ijC ; M ijS g ði; j 2 fr; h; /gÞ are the expansion coefficients of v i ; fi ; rij ; rlij , and M_ ij , respectively. The subscripts C and S, respectively, indicate coefficients of the cosine and sine terms. Similarly, relations between memory variables and components of particle velocity that correspond to Eq. (4) in Cartesian coordinates are described as

^ m Dpl  2Dll @ v^ m @^rlm ðDpl  2Dll Þfm m l @ v rA rrA hA v^ /B ¼ ^r lm þ þ rrA þ Dp r sinh @t @r r @h 2ðDpl  2Dll Þ m Dpl  2Dll m þ v^ rA þ v^ hA coth; ð16Þ r r ^ m Dpl @ v^ m @^rlm ðDpl  2Dll Þfm m l l @ v rA hA  slr hhA ¼ ^rlm v^ /B þ þ hhA þ ðDp  2Dl Þ r sinh @t @r r @h 2ðDpl  Dll Þ m ðDpl  2Dll Þ m þ v^ rA þ v^ hA coth; ð17Þ r r @^rlm ^ m Dpl  2Dll @ v^ m Dpl fm m //A l l @ v rA hA  slr v^ ¼ ^rlm þ þ //A þ ðDp  2Dl Þ r sinh /B @t @r r @h 2ðDpl  Dll Þ m Dpl m þ v^ rA þ v^ hA coth; ð18Þ r r   lm m m ^ @^r 1 @ v^ rA 1 m l @ v hA  slr rhA ¼ ^rlm ð19Þ þ  v^ hA ; rhA þ Dl r @h r @t @r   @^rlm Dll @ v^ m fm m h/B /B v^ hA  v^ m/B coth ; ð20Þ  slr ¼ ^rlm  h/B þ sinh @t r @h  m  lm @^rr/B @ v^ /B fm m 1 m l  slr ð21Þ v^  v^ : ¼ ^rlm  r/B þ Dl r sinh rA r /B @t @r

 slr

lappings (e.g., Virieux, 1986). Seismic waveform calculations by the FDM with standard uniform gridding in spherical coordinates fail near the center of the Earth, because of the following reasons: (1) the extremely small lateral grid spacing near the center perturbs the stability criterion of the FD computations; and (2) there is a singularity of the wave equation in spherical coordinates at the center r ¼ 0. Therefore, our previous work (Toyokuni et al., 2005; Toyokuni and Takenaka, 2006a,b, 2009) retained a doughnut-like computational domain, without a sphere center. For the problem (1), we introduce the multi-domain technique, which uses several computational domains of different grid spacings (see Supplementary information). The multi-domain technique was used in the field of fluid mechanics and then introduced to seismology (e.g., Kessler and Kosloff, 1991; Tessmer et al., 1992; Wang et al., 2001; Wang and Takenaka, 2001; and Wang and Takenaka, 2010 for the PSM, and Falk et al., 1996; Moczo et al., 1996; Aoi and Fujiwara, 1999; Igel and Gudmundsson, 1997; and, Thomas et al., 2000 for the FDM). We constructed a multidomain scheme referring the system of Aoi and Fujiwara (1999) in 3-D Cartesian coordinates, which the upper subdomain to the lower one, with angular grid spacing increased by a factor of three (the 1–3 scheme; see Supplementary information). The FD computations in spherical coordinates require coarser grid spacings in the lateral (h) direction around the center of the sphere, so we should connect several subdomains in the radial (r) direction, to make its lateral grid spacings coarser with depth. Fig. 1 is a schematic of our multi-domain. Denoting the lateral grid spacing of the upper subdomain as Dh, the spacing of lower subdomain becomes 3Dh. We use the FD scheme with second-order accuracy in time and fourth-order accuracy in space (the (2,4) scheme), with a continuous nonuniform grid configuration in the r direction (Pitarka, 1999), so the transmission of variables at the domain boundaries (from lower to upper subdomains) is achieved by interpolation based on fourth-order Lagrange interpolation, to match the spatial order of accuracy with the FDM computation inside each subdomain. For the problem (2), we should consider a scheme to avoid the singularity. Looking back the previous works using the 2.5-D FDM, Thomas et al. (2000) included the Earth center into their calcula-

We have written equations for the cases with Fourier expansion order m – 0, since equations for the m ¼ 0 case can be easily deduced from Eqs. (10)–(21), setting m ¼ 0 and removing A and B. 2.3. Treatment of the Earth center We solve the partial differential equations described above by the velocity-stress, staggered-grid FDM. This calculates three components of the particle velocity vector and six components of the symmetric stress tensor on grid points, allocating different positions for each component of the wavefield, except for several over-

Fig. 1. Schematic of our multi-domain. The number of lateral grid points is trisected from one subdomain to the next deeper one.

49

G. Toyokuni, H. Takenaka / Physics of the Earth and Planetary Interiors 200–201 (2012) 45–55

tions of acoustic wavefields. This time we developed an interpolation technique that enables computations of both P and S waves at the Earth center. Our staggered FD grid has grid points for v h ; v / ; rrr ; rhh ; r// , and rh/ at the center. This forces us to evaluate these wavefield variables exactly on this point, although it is a singularity of wave equations in spherical coordinates. We construct a scheme in which each coefficient for Fourier expansion of the particle velocity and stress at the Earth center is calculated by linear interpolation, using two adjacent values with consideration of symmetry or anti-symmetry of the component about the source axis (see Table 1 of Toyokuni and Takenaka, 2006a):

v^ m ð0; hÞ ¼

v^ m ðDr; hÞ þ v^ m ðDr; h þ pÞ

; 2 r^ m ðDr; hÞ þ r^ m ðDr; h þ pÞ r^ m ð0; hÞ ¼ ; 2

elements in local Cartesian coordinates for source excitation M ij ði; j 2 f1; 2; 3gÞ are related to corresponding elements in global spherical coordinates M kl ðk; l 2 fr; h; /gÞ by equation

0

M rr B @  

M rh Mhh 

0 M33 C B Mh/ A ¼ @  M//  M r/

1

M 13 M 11 

1 M 23 C M 12 A;

ð24Þ

M 22

through assumption that a small grid cell around the source axis in spherical coordinates can be regarded as a rectangular region in Cartesian coordinates.

ð22Þ

3. Results

ð23Þ

In this section, we first check the accuracy and feasibility of our scheme, by comparing the synthetic seismograms with those obtained by the validated method. Then, we show a numerical example of the newly constructed spherical 2.5-D FDM scheme to a realistic Earth model with hemispherical distribution of P-wavespeed and attenuation in the top 300 km of the inner core, using a realistic moment-tensor point source.

where  ¼ 1 for the symmetric, and  ¼ 1 for the anti-symmetric components, respectively. The number of lateral grid points defined at the bottom subdomain is kept even on the center, since the expansion coefficients are functions of angular-dependent unit ba^r and e ^h , which permit putting multiple values of wavesis vectors e fields on this point according to h (as shown in the above equations). 2.4. Wavefield excitation Our FDM scheme can implement an arbitrary moment-tensor point source for excitation of the seismic wavefield. There are two implementations of the seismic source for the velocity–stress, staggered-grid FDM, i.e., implementations on body-force terms (particle velocity implementation, e.g., Frankel, 1993; Graves, 1996), and on moment terms (stress implementation, e.g., Virieux, 1986; Coutant et al., 1995). The two implementations are theoretically equivalent, as supported by the body-force equivalent theorem. As reported in our previous work (Toyokuni and Takenaka, 2006a), we developed a scheme based on the particle velocity implementation. This time, we construct a scheme using the stress implementation, which reduces extent of the source region and is expected to give much accurate waveforms (see Supplementary Figs. S3 and S4). Six moment-tensor elements are implemented so that all elements have the same center of excitation on a grid ^ rr ; r ^ hh ; r ^ // Þ, and for r ^ h/ . The moment-tensor for normal stresses ðr

3.1. Accuracy check To check the accuracy and validity of our FDM scheme, we compare the synthetic seismograms with those obtained by the Direct Solution Method (DSM; Geller and Ohminato, 1994; Takeuchi et al., 1996), which gives exact waveforms for spherically symmetric media. We calculate seismic wave propagation through attenuative Earth structures. First, for simplicity, we superimpose a homogeneous constant Q structure Q P ¼ Q S ¼ 200 on the standard Earth model PREM (Dziewonski and Anderson, 1981) (Accuracy Check # 1; see Supplementary information). When Q P ¼ Q S , strain relaxaSl l tion times for P and S waves become identical (sPl  ¼ s ¼ s ). We compute synthetic seismograms for four situations, with one (L ¼ 1), two (L ¼ 2), five (L ¼ 5), and 10 (L ¼ 10) Zener bodies, to approximate the constant Q model. The stress relaxation times slr are presumed to be reciprocals of the relaxation frequencies covering a frequency range that is logarithmically equidistant. The s-method (Blanch et al., 1995) is then used to optimize strain relaxation times sl for a desired value of quality factor Q ¼ 200 and frequency range ½fa ; fb  (see Supplementary Fig. S5 and Tables S3–S6).

Fig. 2. Synthetic seismograms at Earth surface for the PREM with Q P ¼ Q S ¼ 200. Waveforms of radial component of particle velocity at epicentral distance D ¼ 60 for L ¼ 1; L ¼ 2; L ¼ 5, and L ¼ 10 are displayed. Solid: FDM synthetics; dashed: DSM synthetics. Phases surrounded by black lines in left panel are enlarged in right panel. Observation points for vertical component are located a half grid space below the free surface.

50

G. Toyokuni, H. Takenaka / Physics of the Earth and Planetary Interiors 200–201 (2012) 45–55

Fig. 3. Synthetic seismograms at Earth surface on the / ¼ 0 plane of the PREM with attenuation. Seismic source is a 638-km-deep source with non-zero moment tensor component (a) M 33 ¼ M 0 , (b) M 11 ¼ M 22 ¼ M 0 , (c) M 11 ¼ M 22 ¼ M 0 , (d) M 13 ¼ M 31 ¼ M 0 , (e) M12 ¼ M21 ¼ M0 , and (f) M23 ¼ M32 ¼ M0 . Solid lines: FDM; dashed lines: DSM. Vertical (r) component of particle velocity are shown for (a)–(d), while transverse (/) component are shown for (e) and (f). A lowpass filter (< 1=60 Hz) has been applied to all traces. For the FDM synthetics, observation points for vertical component are located a half grid space below the free surface.

Fig. 2 shows synthetic seismograms calculated by the spherical 2.5-D FDM for the PREM, with Q P ¼ Q S ¼ 200 for the four cases.

This computation used a multi-domain of seven subdomains, with numbers of lateral grids varying from 16 (at the bottom) to 11664

G. Toyokuni, H. Takenaka / Physics of the Earth and Planetary Interiors 200–201 (2012) 45–55

51

Fig. 4. Sequential snapshots of P (red) and S (blue) wave propagation through the PREM, at six time steps after excitation. The region below radius 4000 km has been enlarged to focus on the area around Earth center. The source is a 638-km-deep explosive source. Solid lines are the ICB, CMB, and a border of r ¼ 400 km.

Fig. 5. Source location and mechanism of the 567.8-km-deep Fiji earthquake of 9th March 1994, from the Harvard CMT catalog. Solid and dashed lines indicate a plane including the source and receivers, on which we calculate synthetic seismograms. The solid and dashed lines, respectively, specify eastern ðh P 0Þ and western ðh < 0Þ semicircles about the source axis, as shown in Fig. 6.

(at the top), since the lateral grid spacings triple at every domain interface toward the Earth center. The radial grid spacings Dr and the lateral grid arc length rDh at the Earth surface are 6.0 km and 3.4 km, respectively. The time step is Dt ¼ 0:05 s, and we propagate the signal for 5000 s. To enhance computation accuracy, we have used the effective parameters of density and material parameters in the radial direction, using the GRACE program of Toyokuni and Takenaka (2009). The seismic source is a 635-km-deep dipole source with non-zero moment tensor component M 33 ¼ M0 . The source time function is a phaseless, bell-shaped pulse of width 60 s. Results for L ¼ 1 and L ¼ 2 show faster arrival of succeeding

Fig. 6. Structural model for computation of synthetic seismograms. Thick solid and dashed lines correspond to the same lines in Fig. 5, showing great circle path. We superimpose hemispherical variations of V P and Q P in the top 300 km of the inner core onto the PREM.

phases of the FD synthetics compared with those of the DSM. This is because the desired frequency range for optimization was not adequate for this situation (see Supplementary information). In contrast, results for L ¼ 5 and L ¼ 10 yielded excellent agreement, suggesting that L ¼ 5 is sufficient for an accurate computation. Next, we calculate synthetic seismograms for the PREM with attenuation structure using each moment-tensor element, and compare the results with the DSM synthetics (Accuracy Check # 2; see Supplementary information). We used 10 Zener bodies (L ¼ 10) to approximate the radially varying Q structures. The PREM describes

52

G. Toyokuni, H. Takenaka / Physics of the Earth and Planetary Interiors 200–201 (2012) 45–55

the attenuation structures using quality factor of pure compression Q j and of shear Q l . We obtain Q P and Q S from Q j and Q l using the following relations (p. 187 of Aki and Rechards, 2002):

1 1 ¼ ; QS Ql

  1 1 4V 2 1 1 ; ¼ þ S2  Q P Q j 3V P Q l Q j

ð25Þ ð26Þ

where V P and V S are P- and S-wavespeeds, respectively (see Supplementary Fig. S6). We calculated synthetic seismograms with the same spatial and time grid configuration as the preceding computation. Fig. 3 shows synthetic seismograms at the Earth surface at five epicentral distances on the / ¼ 0 plane, with positive values of h (h P 0) for the moment-tensor element, with nonzero momenttensor components of (a) M 33 , (b) M11 ¼ M 22 , (c) M11 ¼ M22 , (d) M13 ¼ M31 , (e) M12 ¼ M 21 , and (f) M23 ¼ M 32 , respectively. Panels

Fig. 7. Phases related to the inner core at epicentral distances 116°–180 °. Radial component of particle velocity at western ðh < 0Þ semicircle (a) for the PREM; (b) for the inner-core hemispherical heterogeneity; and (c) differential seismograms (b)–(c). The same component at eastern ðh P 0Þ semicircle (d) for the PREM; (e) for the inner-core hemispherical heterogeneity; and (f) differential seismograms (e)–(d). The eastern and western semicircles are indicated in Figs. 5 and 6 using solid and dashed lines, respectively. Travel times of major phases calculated from the TauP Toolkit (Crotwell et al., 1999) for elastic PREM are also marked, using red solid lines. All traces have been low-pass filtered, with corner period of 30 s. Amplitude scale of traces in panels (c) and (f) is twice larger than the other panels.

G. Toyokuni, H. Takenaka / Physics of the Earth and Planetary Interiors 200–201 (2012) 45–55

(a)–(d) show only waveforms of v r , since v / becomes zero everywhere for the axisymmetric elements M 33 and M11 ¼ M22 , and for the double-couple elements M 11 ¼ M 22 and M 13 ¼ M31 ; v / is zero on the / ¼ 0 plane. Owing to a similar reason, panels (e) and (f) only show the v / component of the particle velocity. The fit between the FDM (solid lines) and the DSM (dashed lines) is very good for all source types, phases, and distances, enough to say that these two traces are indistinguishable. We now show sequential snapshots of the seismic wavefield calculated by our method, to demonstrate the feasibility of the multi-domain and interpolation scheme at the Earth center for the PREM. We calculate synthetic seismograms with the same computation parameters described above, but using an explosive source. This source does not have any nodal planes of the wavefield and is therefore useful to test whether waves with large amplitude are accurately propagated through the Earth center. Fig. 4 shows sequential snapshots of P (red) and S (blue) waves, excited by the explosive source and propagating around the Earth center. We have enlarged a region below radius 4000 km, to focus on the center area. According to Wang et al. (2001), P and SV waves are synthesized from the divergence and curl of the wavefield in spherical coordinates, through the following computations.

  1 @ 2  1 @  ðr ur Þ þ ðsin huh Þ; P: 2 r @r r sin h @h   1 @ 1 @ur  ; SV :  ðruh Þ  r @r r @h 

ð27Þ ð28Þ

where ur and uh are r and h components of displacement. At 300 s after excitation, we clearly see in Fig. 4 the incidence of direct P into the fluid outer core. Transmitted core phase PK and converted reflection PcS are propagating toward the center and free surface, respectively. At 450 s, the incidence of succeeding pP into the core produced pPK, and PK traveled through the ICB. At 600 s, the inner core phase PKI have traveled most of the inner core through the center, followed by the converted inner-core shear wave PKJ. At 750 s, PKI has already passed the inner core and formed PKIK, while pPKI and pPKJ are prominent in the inner core. In the frames for 900 and 1050 s, we see in this region a smooth dispersion of P waves from the inner core, as opposed to sustaining reverberation of S waves. These frames show no artificial signals from either the subdomain interfaces or the Earth center, indicating that the current scheme is applicable for realistic Earth models. 3.2. Numerical examples As a numerical example of the present spherical 2.5-D FDM scheme, we calculate the global seismic wavefield for laterally heterogeneous structures in the uppermost inner core with a realistic moment-tensor point source. Many past studies suggest that seismic heterogeneities exist in the top of the inner core (e.g., Kaneshima, 1996; Tanaka and Hamaguchi, 1997). Among others, we deal with the hemispherical variation of P-wavespeed and Q P in this region for our numerical example. Niu and Wen (2001) have analyzed differential PKiKP–PKIKP travel times. They found that within the top 100 km inner core, residuals are systematically larger, by about 0.8 s, for waves sampled in the eastern hemisphere (40°E to 180°E) compared with those sampled in the western hemisphere (180°W to 40°E). Wen and Niu (2002) did an extensive search for these phases, revealing that the faster eastern hemisphere has lower Q than the slower western hemisphere. They showed that the top of eastern and western hemispheres have average Q P values of 250 and 600, respectively. Cao and Romanowicz (2004) also reported large differences in Q P in the uppermost inner core, with Q P  160 in the eastern, and Q P  335 in the western hemisphere. Yu and Wen (2006) did a joint analysis,

53

with observations of PKiKP/PKIKP in the distance range 120°– 141°, and PKP(BC)/PKIKP in the distance range 146°–160°, along equatorial paths. They were able to track the hemispherical differences in Q down to at least 200 km below the ICB. To the contrary, Li and Cormier (2002) analyzed PKIKP data at longer distances, concluding that there was no evidence for such hemispherical variations in the deeper inner core. For a numerical example of our scheme, we construct an Earth model with hemispherical inner core heterogeneity superimposed on the isotropic PREM, without the ocean layer (Dziewonski and Anderson, 1981). We model the 567.8-km-deep Fiji earthquake of 9 March 1994 that occurred at (178:11 W, 17:69 S), which had a magnitude of Mw ¼ 7:6. We used the source mechanism from the Harvard catalog (see Supplementary Table S7). Figs. 5 and 6 indicate the cross-section for FDM computation. It is a nearly equatorial plane including the seismic source, and a boundary between the eastern and western inner-core hemispheres is selected with axes h ¼ 0 and h ¼ p. In the top 300 km of the inner core, the eastern hemisphere ð1:89 E to 178:11 W through 180 E) has a +20% faster V P value and Q P ¼ 300, whereas the western hemisphere ð178:11 W to 1:89 E through 0 E) has a 20% slower V P value and Q P ¼ 600 (Fig. 6). We construct a multi-domain composed of seven subdomains, with numbers of lateral grid points varying from 32 (at the bottom) to 23328 (at the top) (see Supplementary Table S8). At the free surface, radial and lateral grid spacings are about 3.3 km and 1.7 km, respectively. The number of spatial grid points on a cross-section of the Earth is about 1  107 . The time step is Dt ¼ 0:025 s, and we propagate the signal for 5000 s. Fig. 7 shows synthetic seismograms calculated for the PREM (panels (a) and (d)) and for the model with a laterally heterogeneous inner core (panels (b) and (e)), as well as differential seismograms obtained from subtraction of the PREM synthetics from the synthetics for the laterally heterogeneous inner core (panels (c) and (f)). Panels (a), (b), and (c) are for h < 0, i.e., the western semicircle indicated by dashed lines in Figs. 5 and 6. The other three panels, (d), (e), and (f), are for h P 0, i.e., the eastern semicircle represented by solid lines in these figures. These panels show arrival times of major inner-core phases calculated by the elastic PREM, using the TauP toolkit (Crotwell et al., 1999). Panels (c) and (f) clearly show that most seismic phases affected by the lateral heterogeneity of the inner core are inner-core phases, as expected. Arrivals of S-related phases, e.g., SKIKS and sSKIKS, are more delayed from expected arrival times with the elastic model, in comparison with P-related phases. This is because the PREM imposes larger attenuation on S waves. We may confirm that the I-related phases, such as PKIKP and pPKIKP, arrive faster (by a maximum of 20 s) in the western semicircle than in the eastern semicircle, because of the 20% V P heterogeneity in the top 300 km of the PREM inner core. This result also shows interesting characteristics of core phases. Namely, in the eastern semicircle, the differential amplitudes of the I-related phases, especially pPKIKP, are apparently enhanced by a concentration of wavefields, whereas those in the western semicircle tend to be dispersed.

4. Discussion and conclusions This paper had provided a method for realistic global seismic waveform modeling using an arbitrary, laterally heterogeneous, attenuative Earth model including the Earth center, and an arbitrary moment-tensor point source. The method is based on spherical 2.5-D modeling that solves the 3-D viscoelastodynamic equation in spherical coordinates, on a cross-section along a great circle of the Earth including the source and receivers. It requires computation resources as small as 2-D modeling since it assumes that the majority of seismic energy is transported on this plane.

54

G. Toyokuni, H. Takenaka / Physics of the Earth and Planetary Interiors 200–201 (2012) 45–55

However, the method can correctly account for the 3-D pulse shape and geometrical spreading effects, in contrast to 2-D modeling, enabling direct comparison of synthetic seismograms with observations. In this instance, we have applied a scheme using memory variables to incorporate anelastic attenuation into our timedomain, velocity–stress, staggered-grid FD computations. The memory-variable schemes have been improved in Cartesian coordinates for decades, but this is the first application to the FD computations in spherical coordinates. We have confirmed that the scheme has sufficient accuracy for the current situation, by comparing the results with the DSM synthetics at the free surface of the Earth. Further precise approaches have been recently proposed, such as material-independent memory variables (Kristek and Moczo, 2003), which only allow for spatial coarse sampling of these variables in the presence of a material discontinuity. Implementation of such up-to-date schemes will be the subject of future work. We have also solved problems related to the Earth center, through the multi-domain technique, which has wider lateral grid spacings in the deep part of the Earth, and through interpolation of wavefield variables exactly on the center, considering symmetry and anti-symmetry of wave distribution about the source axis with respect to a specific moment-tensor element. The current scheme has applied to modeling global seismic waveforms excited from a point source with a realistic mechanism, which propagate through a realistic Earth with hemispherical variation V P and Q P in the top 300 km of the inner core, using a source time function of duration 30 s. The calculation requires five FD computations of the Fourier expansion coefficients for a moment-tensor element. For a total duration of 5000 s after excitation, each FD computation required 5.6 GB of memory in a single precision calculation and a computation time of 4.2 days, using eight CPUs with IBM POWER6 architecture (4.7 GHz clock speed). This is one thousandth or even one ten-thousandth of the resources required for 3-D simulation. Making the best use of this advantage, Toyokuni et al. (2012) calculated global SH wave propagation up to the period of 4 s, which is useful for investigating short-period core-reflected S waves. Recently, Wookey and Helffrich (2008) proposed that the inner core Q S might be higher than previously estimated. As a near-future work, comparison of high-frequency synthetics of the full wavefield for the different inner core Q models could be suitable for our approach. Although it requires much computation time compared with quasi-analytic method such as the reflectivity method, our scheme can calculate complete wavefield including all reflected, transmitted, and converted body waves, as well as surface waves. Furthermore, it simultaneously provides wavefield values at all spatial grid points spread over the cross section, so we can intuitively grasp the nature of wave propagation through, for example, snapshots. These resulting wave traces can even be directly compared with observations. For example, synthetic seismograms calculated by our method showed excellent agreement with three-component observed seismograms of the November 9, 2009, 603.9 km deep Fiji earthquake (Toyokuni et al., under review). In global seismology, recent increases in data and resolution of global and regional structures have accelerated theoretical works toward the use of complete 3-D modeling, entailing vast simulations of the Earth’s inner structure. Our scheme represents another route toward this goal. A proven and efficient technique could easily serve as a workhorse in the iterative calculation of synthetics for waveform inversion, or for adjusting an initial model of 3-D computation. Acknowledgments We are grateful to Prof. George Helffrich and two anonymous reviewers for kind and constructive review. Most figures were

generated using the Generic Mapping Tools freeware package (Wessel and Smith, 1998). This research was partly supported by JSPS.KAKENHI(189721, 216753). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.pepi.2012.03.010. References Aki, K., Richards, P.G., 2002. Quantitative Seismology, second ed. University Science Books, Sausalito, CA, 700 pp. Alterman, Z.S., Aboudi, J., Karal, F.C., 1970. Pulse propagation in a laterally heterogeneous solid elastic sphere. Geophys. J. R. Astr. Soc. 21 (3), 243–260. Aoi, S., Fujiwara, H., 1999. 3D finite-difference method using discontinuous grids. Bull. Seism. Soc. Am. 89 (4), 918–930. Blanch, J.O., Robertsson, J.O.A., Symes, W.W., 1995. Modeling of a constant Q: methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique. Geophysics 60 (1), 176–184. Cao, A., Romanowicz, B., 2004. Hemispherical transition of seismic attenuation at the top of the earth’s inner core. Earth Planet. Sci. Lett. 228 (3–4), 243–253. http://dx.doi.org/10.1016/j.epsl.2004.09.032. Carcione, J.M., Kosloff, D., Kosloff, R., 1988a. Wave propagation simulation in a linear viscoelastic medium. Geophys. J. 95 (3), 597–611. Carcione, J.M., Kosloff, D., Kosloff, R., 1988b. Viscoacoustic wave propagation simulation in the earth. Geophysics 53 (6), 769–777. Chaljub, E., Tarantola, A., 1997. Sensitivity of SS precursors to topography on the upper-mantle 660-km discontinuity. Geophys. Res. Lett. 24 (21), 2613–2616. Chaljub, E., Capdeville, Y., Vilotte, J.-P., 2003. Solving elastodynamics in a fluid-solid heterogeneous sphere: a parallel spectral element approximation on nonconforming grids. J. Comp. Phys. 187 (2), 457–491. http://dx.doi.org/10.1016/ S0021-9991(03)00119-0. Coutant, O., Virieux, J., Zollo, A., 1995. Numerical source implementation in a 2D finite difference scheme for wave propagation. Bull. Seism. Soc. Am. 85 (5), 1507–1512. Crotwell, H.P., Owens, T.J., Ritsema, J., 1999. The TauP toolkit: flexible seismic travel-time and ray-path utilities. Seism. Res. Lett. 70 (2), 154–160. http:// dx.doi.org/10.1785/gssrl.70.2.154. Dziewonski, A.M., Anderson, D.L., 1981. Preliminary reference Earth model. Phys. Earth Planet. Int. 25 (4), 297–356. Emmerich, H., Korn, M., 1987. Incorporation of attenuation into time-domain computations of seismic wave fields. Geophysics 52 (9), 1252–1264. Falk, J., Tessmer, E., Gajewski, D., 1996. Tube wave modeling by the finite-difference method with varying grid spacing. PAGEOPH 148 (1–2), 77–93. Frankel, A., 1993. Three-dimensional simulations of ground motions in the San Bernardino Valley California for hypothetical earthquakes on the San Andreas fault. Bull. Seism. Soc. Am. 83 (4), 1020–1041. Geller, R.J., Hara, T., 1993. Two efficient algorithms for iterative linearized inversion of seismic waveform data. Geophys. J. Int. 115 (3), 699–710. http://dx.doi.org/ 10.1111/j.1365-246X.1993.tb01488.x. Geller, R.J., Ohminato, T., 1994. Computation of synthetic seismograms and their partial derivatives for heterogeneous media with arbitrary natural boundary conditions using the Direct Solution Method. Geophys. J. Int. 116 (2), 421–446. Graves, R.W., 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences. Bull. Seism. Soc. Am. 86 (4), 1091–1106. Igel, H., Gudmundsson, O., 1997. Frequency-dependent effects on travel times and waveforms of long-period S and SS waves. Phys. Earth Planet. Int. 104 (1–3), 229–246. Igel, H., Weber, M., 1995. SH-wave propagation in the whole mantle using highorder finite differences. Geophys. Res. Lett. 22 (6), 731–734. Igel, H., Weber, M., 1996. P-SV wave propagation in the Earth’s mantle using finite differences: application to heterogeneous lowermost mantle structure. Geophys. Res. Lett. 23 (5), 415–418. JafarGandomi, A., Takenaka, H., 2007. Efficient FDTD algorithm for plane-wave simulation for vertically heterogeneous attenuative media. Geophysics 72 (4), H43–H53. Jahnke, G., Thorne, M.S., Cochard, A., Igel, H., 2008. Global SH-wave propagation using a parallel axisymmetric spherical finite-difference scheme: application to whole mantle scattering. Geophys. J. Int. 173 (3), 815–826. http://dx.doi.org/ 10.1111/j.1365-246X.2008.03744.x. Kaneshima, S., 1996. Mapping heterogeneity of the uppermost inner core using two pairs of core phases. Geophys. Res. Lett. 23 (22), 3075–3078. Kessler, D., Kosloff, D., 1991. Elastic wave propagation using cylindrical coordinates. Geophysics 56 (12), 2080–2089. Komatitsch, D., Tromp, J., 2002a. Spectral-element simulations of global seismic wave propagation—I. Validation. Geophys. J. Int. 149 (2), 390–412. Komatitsch, D., Tromp, J., 2002b. Spectral-element simulations of global seismic wave propagation—II. Three-dimensional models oceans rotation and selfgravitation. Geophys. J. Int. 150 (1), 303–318. http://dx.doi.org/10.1046/j.1365246X.2002.01716.x. Komatitsch, D., Michéa, D., Erlebacher, G., 2009. Porting a high-order finite-element earthquake modeling application to NVIDIA graphics cards using CUDA. J.

G. Toyokuni, H. Takenaka / Physics of the Earth and Planetary Interiors 200–201 (2012) 45–55 Parallel Distrib. Comput. 69 (5), 451–460. http://dx.doi.org/10.1016/j.jpdc. 2009.01.006. Kristek, J., Moczo, P., 2003. Seismic-wave propagation in viscoelastic media with material discontinuities a 3D fourth-order staggered-grid finite-difference modeling. Bull. Seism. Soc. Am. 93 (5), 2273–2280. http://dx.doi.org/10.1785/ 0120030023. Li, X., Cormier, V.F., 2002. Frequency-dependent seismic attenuation in the inner core 1. A viscoelastic interpretation. J. Geophys. Res. 107 (B12), 2361. http:// dx.doi.org/10.1029/2002JB001795. Liu, H.-P., Anderson, D.L., Kanamori, H., 1976. Velocity dispersion due to anelasticity; implications for seismology and mantle composition. Geophys. J. R. Astr. Soc. 47 (1), 41–58. Moczo, P., Kristek, J., 2005. On the rheological models used for time-domain methods of seismic wave propagation. Geophys. Res. Lett. 32, L01306. http:// dx.doi.org/10.1029/2004GL021598. Moczo, P., Kristek, J., Halada, L., 2004. The Finite-Difference Method for Seismologists. An Introduction. Comenius University, Bratislava, 158 pp., ISSN: 80-223-2000-5. Moczo, P., Labák, P., Kristek, J., Hron, F., 1996. Amplification and differential motion due to an antiplane 2D resonance in the sediment valleys embedded in a layer over the half-space. Bull. Seism. Soc. Am. 86 (5), 1434–1446. Nissen-Meyer, T., Fournier, A., Dahlen, F.A., 2007. A two-dimensional spectralelement method for computing spherical-earth seismograms—I. Momenttensor source. Geophys. J. Int. 168 (3), 1067–1092. http://dx.doi.org/10.1111/ j.1365-246X.2006.03121.x. Nissen-Meyer, T., Fournier, A., Dahlen, F.A., 2008. A 2-D spectral-element method for computing spherical-earth seismograms—II. Waves in solid-fluid media. Geophys. J. Int. 174 (3), 873–888. http://dx.doi.org/10.1111/j.1365-246X.2008. 03813.x. Niu, F., Wen, L., 2001. Hemispherical variations in seismic velocity at the top of the Earth’s inner core. Nature 410, 1081–1084. http://dx.doi.org/10.1038/ 35074073. Okamoto, T., Takenaka, H., Nakamura, T., Aoki, T., 2010. Accelerating large-scale simulation of seismic wave propagation by multi-GPUs and three-dimensional domain decomposition. Earth Planets Space 62 (12), 939–942. http://dx.doi.org/ 10.5047/eps.2010.11.009. Pitarka, A., 1999. 3D elastic finite-difference modeling of seismic motion using staggered grids with nonuniform spacing. Bull. Seism. Soc. Am. 89 (1), 54–68. Robertsson, J.O.A., Blanch, J.O., Symes, W.W., 1994. Viscoelastic finite-difference modeling. Geophysics 59 (9), 1444–1456. http://dx.doi.org/10.1190/1.1443701. Takenaka, H., Furumura, T., Fujiwara, H., 1998. Recent developments in numerical methods for ground motion simulation. In: Irikura, K., Kudo, K., Okada, H., Sasatani, T. (Eds.), The Effects of Surface Geology on Seismic Motion, vol. 1. Balkema, Rotterdam (ISBN: 90-5809-030-2), pp. 91–101. Takeuchi, N., Geller, R.J., Cummins, P.R., 1996. Highly accurate P-SV complete synthetic seismograms using modified DSM operators. Geophys. Res. Lett. 23 (10), 1175–1178. Tanaka, S., Hamaguchi, H., 1997. Degree one heterogeneity and hemispherical variation of anisotropy in the inner core from PKP(BC)–PKP(DF) times. J. Geophys. Res. 102 (B2), 2925–2938. Tessmer, E., Kessler, D., Kosloff, D., Behle, A., 1992. Multi-domain Chebyshev– Fourier method for the solution of the equations of motion of dynamic elasticity. J. Comput. Phys. 100 (2), 355–363. http://dx.doi.org/10.1016/00219991(92)90241-P. Thomas, Ch., Igel, H., Weber, M., Scherbaum, F., 2000. Acoustic simulation of P-wave propagation in a heterogeneous spherical earth: numerical method and

55

application to precursor waves to PKPdf. Geophys. J. Int. 141 (2), 307–320. http://dx.doi.org/10.1046/j.1365-246x.2000.00079.x. Thorne, M.S., Lay, T., Garnero, E.J., Jahnke, G., Igel, H., 2007. Seismic imaging of the laterally varying D00 region beneath the cocos plate. Geophys. J. Int. 170 (2), 635– 648. http://dx.doi.org/10.1111/j.1365-246X.2006.03279.x. Toyokuni, G., Takenaka, H., 2006a. FDM computation of seismic wavefield for an axisymmetric earth with a moment tensor point source. Earth Planets Space 58 (8), e29–e32. Toyokuni, G., Takenaka, H., 2006b. Efficient method to model seismic wave propagation through deep inside the earth: quasi-spherical approach. Chikyu Monthly 28 (9), 607–611 (in Japanese). Toyokuni, G., Takenaka, H., 2009. ACE–A FORTRAN subroutine for analytical computation of effective grid parameters for finite-difference seismic waveform modeling with standard Earth models. Comput. Geosci. 35 (3), 635–643. http://dx.doi.org/10.1016/j.cageo.2008.05.005. Toyokuni, G., Takenaka, H., Kanao, M., 2012. Quasi-axisymmetric finite-difference method for realistic modeling of regional and global seismic wavefield – review and application. In: M. Kanao (Ed.), Seismic Waves, Research and Analysis, InTech (ISBN: 978-953-307-944-8), pp. 85–112. Toyokuni, G., Takenaka, H., Wang, Y., Kennett, B.L.N., 2005. Quasi-spherical approach for seismic wave modeling in a 2-D slice of a global Earth model with lateral heterogeneity. Geophys. Res. Lett. 32, L09305. http://dx.doi.org/ 10.1029/2004GL022180. Toyokuni, G., Takenaka, H., Kanao, M., Wiens, D., Nyblade, A., under review. Comparison of global synthetic seismograms calculated by the spherical 2.5-D finite-difference method with observed long-period waveforms including data from intra-Antarctic region. Polar Sci.. Tsuboi, S., Komatitsch, D., Ji, C., Tromp, J., 2008. Computations of global seismic wave propagation in three dimensional earth model. In: Labarta, J., Joe, K., Sato, T. (Eds.), Proceedings of the Sixth International Symposium on High-performance Computing and First International Conference on Advanced Low Power Systems, vol. 4759. Springer, Berlin/Heidelberg (ISBN: 978-3-540-77703-8), pp. 434–443. Virieux, J., 1986. P-SV wave propagation in heterogeneous media: velocity–stress finite-difference method. Geophysics 51 (4), 889–901. Wang, Y., Takenaka, H., 2001. A multidomain approach of the Fourier pseudospectral method using discontinuous grid for elastic wave modeling. Earth Planets Space 53 (3), 149–158. Wang, Y., Takenaka, H., 2010. A scheme to treat the singularity in global seismic wavefield simulation using pseudospectral method with staggered grids. Earthquake Sci. 23 (2), 121–127. http://dx.doi.org/10.1007/s11589-010-0001-x. Wang, Y., Takenaka, H., Furumura, T., 2001. Modelling seismic wave propagation in a two-dimensional cylindrical whole-earth model using the pseudospectral method. Geophys. J. Int. 145 (3), 689–708. http://dx.doi.org/10.1046/j.1365246x.2001.01413.x. Wen, L., Niu, F., 2002. Seismic velocity and attenuation structures in the top of the Earth’s inner core. J. Geophys. Res. 107 (B11), 2273. http://dx.doi.org/10.1029/ 2001JB000170. Wessel, P., Smith, W.H.F., 1998. New, improved version of generic mapping tools released. Eos Trans. AGU 79 (47), 579. http://dx.doi.org/10.1029/98EO00426. Wookey, J., Helffrich, G., 2008. Inner-core shear-wave anisotropy and texture from an observation of PKJKP waves. Nature 454, 873–877. http://dx.doi.org/ 10.1038/nature07131. Yu, W., Wen, L., 2006. Seismic velocity and attenuation structures in the top 400 km of the Earth’s inner core along equatorial paths. J. Geophys. Res. 111, B07308. http://dx.doi.org/10.1029/2005JB003995.