Quasi-Axisymmetric Finite-Difference Method for Realistic Modeling of Regional and Global Seismic Wavefield — Review and Application —

In this chapter, we describe recent developments of forward-modeling techniques for accurate and efficient computation of the realistic seismic wavefield. Our knowledge on the Earth’s interior has been enhanced by mutual progress in observation and numerical methods. Since the first time-recording seismograph was built in Italy in 1875 (Shearer, 1999), the recorded seismic dataset has been growing at an almost exponential rate. Such a massive amount of seismic waveform data should be interpreted with consideration of the seismic source mechanism and Earth’s inner structure, which explain each crest or trough in observed waveform traces. This interpretation can be achieved by forward modeling of seismic waveforms. In addition, recent progress in computation capacity has enabled investigation of the Earth’s inner structure via waveform inversion, an inverse problem minimizing the difference between observed and synthetic seismograms. This method requires iterative computations of synthetic seismograms for each structural model renewal in theminimization process, so we need a forward modeling technique that produces accurate waveforms with small computation time and memory.


Introduction
In this chapter, we describe recent developments of forward-modeling techniques for accurate and efficient computation of the realistic seismic wavefield.Our knowledge on the Earth's interior has been enhanced by mutual progress in observation and numerical methods.Since the first time-recording seismograph was built in Italy in 1875 (Shearer, 1999), the recorded seismic dataset has been growing at an almost exponential rate.Such a massive amount of seismic waveform data should be interpreted with consideration of the seismic source mechanism and Earth's inner structure, which explain each crest or trough in observed waveform traces.This interpretation can be achieved by forward modeling of seismic waveforms.In addition, recent progress in computation capacity has enabled investigation of the Earth's inner structure via waveform inversion, an inverse problem minimizing the difference between observed and synthetic seismograms.This method requires iterative computations of synthetic seismograms for each structural model renewal in the minimization process, so we need a forward modeling technique that produces accurate waveforms with small computation time and memory.
Writing mathematically, forward modeling (forward problem, modelization problem, or simulation problem) predicts error-free values of observable parameters d corresponding to a given model m, i.e., this theoretical prediction can be denoted where d = g(m) is a short notation for a set of equations The operator g(•) is called the forward operator, which expresses our mathematical model of the physical system under study (Tarantola, 2005).The forward modeling of seismic waveforms is therefore a theoretical method that applies a set of theoretical equations to determine what given seismographs would measure with respect to a preset combination of source and structure.Basically, the forward modeling of seismic waves solves the elastodynamic equation for a given source mechanism and structural model, including a set of density and elastic parameters.
advantage of the staggered grid is its robustness for structures with high contrast of Poisson's ratio.
The elastodynamic equation (equations of motion and constitutive laws), together with the initial and boundary conditions, completely describe a problem of seismic wave propagation.
If we keep the equations of motion separate from the constitutive laws, we can speak of the displacement-stress formulation.If we use particle velocity in the equations of motion, keep the constitutive laws, and add the definition of particle velocity, we obtain the displacement-velocity-stress formulation.If we apply a time derivative to the constitutive laws instead of adding the definition of particle velocity, then we have the velocity-stress formulation.If we eliminate the stress-tensor components by substituting the constitutive laws into the equations of motion, we get the displacement formulation (e.g., Moczo et al., 2007).In this chapter, we mainly deal with the velocity-stress formulation.The formulation for 3-D computations for general elastic media are described with mathematical expressions in Cartesian coordinates (x 1 , x 2 , x 3 ) as follows (i, j, k, l ∈{1, 2, 3}): where t is time, ρ(x); x =( x 1 , x 2 , x 3 ) is the density; c ijkl (x) is the component of the elastic tensor.In addition, v i (x, t), f i (x, t), σ ij (x, t), ǫ ij (x, t),a n d Ṁij (x, t) are components of the particle velocity vector, body force vector, stress tensor, strain tensor, and first-order time derivative of the moment tensor, respectively.We have used the summation convention over repeated suffixes.
In general orthogonal curvilinear coordinates (c 1 , c 2 , c 3 ), the elastodynamic equation corresponds to Eqs. (2), and (3) can be given as follows, together with the strain-velocity relation without using the summation convention (Aki & Richards, 2002): where h i ; i ∈{1, 2, 3} are scaling functions peculiar to the coordinate system, and all variables related to the medium and wavefields are the same as in the Cartesian case, except that the x dependence should be replaced with c =(c 1 , c 2 , c 3 ) dependence.We consider special cases of these equations in cylindrical and spherical coordinates in the following sections.
With simulation of seismic wave motion by the domain methods, two "dimensions", i.e., the spatial dimension (heterogeneity) of a medium and the dimensionality of wavefields, become important.The heterogeneity of a medium is defined as the number of independent 87 Quasi-Axisymmetric Finite-Difference Method for Realistic Modeling of Regional and Global Seismic Wavefield -Review and Applicationwww.intechopen.comvariables considering material parameters as spatial functions, whereas the dimensionality of wavefields is defined as the number of spatial variables in all independent variables of wavefields.We should choose these two dimensions for a compromise between accuracy and efficiency of waveform computation.Takenaka (1993;1995) developed a notation describing the combination as (m, n) dimension, where m is the dimensionality of wavefields and n is the dimension of a medium.The (3, 3) dimensional modeling, called 3-D modeling, calculates 3-D seismic wavefields in a 3-D structural model.The 3-D modeling provides accurate results, since it can treat the most realistic situation.However, full 3-D modeling up to a realistic high frequency is still computationally intensive and costly, even on parallel hardware.On the other hand, (2, 2) dimensional modeling (2-D modeling) calculates 2-D wavefields in a 2-D structural model, which requires relatively small computational resources compared to 3-D modeling.Because of the severe computational requirement, waveform computation by domain methods had long been restricted to 2-D modeling, although the wave behavior of out-of-plane motion is underestimated.Moreover, 2-D modeling cannot correctly model geometrical spreading effects and the pulse shape in 3-D.
In the 1990s, an alternative (3, 2) dimensional modeling called 2.5-D modeling, which calculates 3-D wavefields in a medium varying only in two dimensions, was introduced in seismology.It first assumes the structure to be invariant in one direction, and then applies a spatial Fourier transform to the 3-D wave equation in this direction.The resulting equations in a mixed coordinate-wavenumber domain consist of independent sets of 2-D equations for each wavenumber, such that numerical computations of these equations followed by an inverse Fourier transform over wavenumber generate 3-D synthetic seismograms.One can thus correctly model the 3-D geometrical spreading effects and pulse shape for all phases, and it makes possible a direct comparison between real and synthetic waveform data.Nevertheless, associated with computations for all discrete wavenumbers, 2.5-D modeling requires a long computation time, comparable to 3-D modeling, although it has a memory requirement only slightly greater than 2-D modeling.A more economical technique for modeling 3-D seismic wavefields is to approximate the structural model as rotationally symmetric along the vertical axis, include a seismic source, and then solve the wave equation in cylindrical or spherical coordinates.This method,

88
Seismic Waves, Research and Analysis called axisymmetric modeling, can correctly model 3-D geometrical spreading effects and pulse shape, with computation time and memory comparable to 2-D modeling.In the following subsections, we write the elastodynamic equations for axisymmetric modeling in both cylindrical coordinates (r, φ, z) and spherical coordinates (r, θ, φ), as shown in Figure 1, comparing them with the 3-D equations.

For axisymmetric modeling
Axisymmetric modeling in cylindrical coordinates assumes the structure to be axisymmetric with respect to the axis r = 0.For cases with axisymmetric seismic sources, such as explosive and torque sources (SH-wave source), the 3-D seismic wavefield is completely separated into in-plane motion in the r-z plane (P-SV waves) and anti-plane motion (SH waves).In this situation, terms including φ derivatives can be neglected in Eqs. ( 8)-( 16), which gives, for example, the P-SV elastodynamic equation for axisymmetric modeling in cylindrical coordinates from axisymmetric sources as

For axisymmetric modeling
Axisymmetric modeling in spherical coordinates assumes the structure to be axisymmetric about the axis θ = 0.As we have already seen in Section 3.1.2,the in-plane P-SV motion in the r-θ plane and the anti-plane SH motion can completely been separated for axisymmetric sources.Consequently, the P-SV elastodynamic equation for axisymmetric modeling in spherical coordinates from axisymmetric sources becomes Similarly, the SH elastodynamic equation in spherical coordinates for axisymmetric sources becomes the following.
This decoupling between P-SV and SH waves only holds for axisymmetric sources.Toyokuni & Takenaka (2006a) implemented arbitrary moment-tensor point sources, including shear dislocation sources, into the axisymmetric computation in spherical coordinates, using the Fourier transform of all field variables in the φ direction, which can be written as where a is a variable that can be replaced by any component of the particle velocity vector, body force vector, stress tensor, and moment tensor; m is the expansion order and â0 , âm C , âm S are expansion coefficients.Subscripts C and S have been added to indicate coefficients for cosine and sine terms, respectively.It is sufficient to take the expansion order up to m = 2 with consideration of radiation patterns of moment tensor sources.Substitution of Eq. ( 42) into the 3-D elastodynamic equation in spherical coordinates Eqs. ( 24)-( 32), followed by rearrangement, gives five closed systems of the partial differential equations of the expansion coefficients.The equations for m = 0 have the same form as Eqs.( 33)-( 41), whereas those for m = 1, 2 become the following: where A, B ∈{ C, S}.Note that P-SV and SH waves are coupled via coupling terms in these equations, including m.Through the Fourier expansion in the φ direction, an arbitrary moment tensor source is decomposed into five moment tensor elements: (1) Axisymmetric excitation for m = 0, and four double couple excitations for m = 1, 2 that are classified using the combinations of three parameters {m, A, B} as (2) {1, C, S}, (3) {1, S, C}, (4) {2, C, S}, and (5) {2, S, C}.The elements (2) and (3) correspond to purely vertical dip-slip excitations shifted π/2 in the φ direction for each other, whereas (4) and ( 5) are purely strike-slip excitations shifted π/4 in the φ direction for each other, as shown in Figure 2. Computations of expansion coefficients using, for example, the FDM with respect to each moment tensor element via Eqs.( 33)-( 41) or Eqs. ( 43)-( 51), followed by substitution of the results into Eq.( 42), enables attainment of global elastic response by an arbitrary moment tensor source for the axisymmetric structural model.This requires only five times the computational resources of computations for purely axisymmetric sources.

Review of axisymmetric modeling
Because of the light computational requirement and correct treatment of 3-D seismic wavefields, axisymmetric modeling has frequently been used by researchers.Here, we briefly summarize their works.

92
Seismic Waves, Research and Analysis Axisymmetric modeling in cylindrical coordinates has often been used to efficiently calculate realistic 3-D seismic wavefields, especially for target areas of seismic exploration.For flat-lying media, the solution on a r-z cross section with a source and receivers will be correct for full 3-D modeling for a point source.When using Cartesian 2-D modeling for the same target, the seismic source becomes a line in 3-D (point in 2-D) along a direction of structural invariance.This causes fatal errors on waveforms and makes it impossible for direct comparison between real and synthetic data, even when the real data are converted to 2-D.However, in axisymmetric modeling in cylindrical coordinates, any lateral variations on the r-z plane become physically unrealistic rings in 3-D, except in very special cases.Alterman & Karal (1968) introduced the technique for FDM computation of seismic waves in elastic layered half-space with a buried point compressional source.They applied the scheme to various investigations, e.g., of the effect of different mesh sizes on synthetic seismograms, development of Rayleigh waves on the surface, change of Rayleigh waves with depth and pulse width, and so on.Details of their FDM scheme are also in Alterman & Loewenthal (1972).Stephen (1983) adopted the cylindrical approach to compare the FDM and reflectivity synthetic seismograms for a compressional point source, using laterally homogeneous seafloor models with step and ramp discontinuities between liquid and solid, showing the two methods to be in excellent agreement.Stephen (1988) expanded the work of Stephen (1983), testing various FDM formulations to determine which ones produce acceptable solutions for seafloor problems.They used models with horizontal liquid-solid interfaces, and those with rough shape.Igel et al. (1996) performed waveform inversion of marine reflection seismograms to determine P impedance and Poisson's ratio structures in the Gulf of Mexico, through iterative calculation of synthetic seismograms by axisymmetric modeling in cylindrical coordinates.
Axisymmetric modeling in spherical coordinates is a powerful technique to obtain the realistic 3-D global seismic wavefield.Therefore, it has long been used, in spite of the restriction of structural models in rotational symmetry with respect to the axis through the seismic source.Alterman et al. (1970) pioneered the application of this method to the FDM computations of elastic wavefield radiated by an impulsive point source, for radially and laterally heterogeneous, purely mathematical sphere models.The first application of this approach to the Earth model was the work of Igel & Weber (1995).They simulated SH-wave propagation in frequency bands up to 0.1 Hz in the whole mantle model from the PREM (Dziewonski & Anderson, 1981), using displacement-stress FD schemes with an eight-point 93 Quasi-Axisymmetric Finite-Difference Method for Realistic Modeling of Regional and Global Seismic Wavefield -Review and Applicationwww.intechopen.comoperator in space.Igel & Weber (1996) extended this high-order FDM scheme to P-SV waves.They investigated the effects of heterogeneities in the D" layer on long-period P-waves, with a dominant period of 15 s excited by an explosive source.Their computational domain was restricted to a region with angular range of 105 • and maximum depth 4600 km, because of high computational requirements and perturbation of the FD stability criterion near the Earth center.Nevertheless, they successfully examined three lower mantle models, in addition to the isotropic PREM.Chaljub & Tarantola (1997) also used the displacement-stress FDM scheme for SH waves, to test sensitivity of SS precursors to the presence of topography on the 660-km discontinuity.They adopted models with topography on the discontinuity, as well as with a penetrating slab toward it at various scales, to examine its apparent depth deduced by bottomside reflection S660S.Igel & Gudmundsson (1997) applied a multi-domain, i.e., the FDM grid configuration with vertically-varying lateral grid spacing, to the SH algorithm.This was done to investigate frequency-dependent effects on S and SS waveforms and travel times through random upper-mantle models with pre-assumed spectral properties.Thomas et al. (2000) solved the acoustic wave equation by axisymmetric modeling in spherical coordinates, using the multi-domain including the Earth center.They used the scheme to study the influence of velocity contrasts, location, and orientation of various scatterers imposed near the CMB on precursors to PKPdf.Although axisymmetric modeling itself can treat an arbitrary moment-tensor point source, all works listed above concentrated on using axisymmetric sources, such as explosive and torque sources.Toyokuni & Takenaka (2006a) therefore developed a scheme to implement a non-axisymmetric source in the FDM scheme based on axisymmetric modeling in spherical coordinates, using the Fourier expansion of wavefield variables in the φ direction, as in Section 3.2.2.As a numerical example, they simulated which seismic phases can be related to a stagnant slab located far from a point source, with the mechanism of the 1994 deep Bolivia earthquake.Jahnke et al. (2008) extended the SH scheme of Igel & 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.

Quasi-axisymmetric modeling
As stated in the previous section, axisymmetric modeling remains a powerful tool to obtain the 3-D seismic wavefield, because its economical calculation focuses only on a cross section, including the source and receivers.Especially in global modeling, axisymmetric modeling in spherical coordinates is the best way for iterative computation of synthetic seismograms for inverting data to image the Earth's inner structures by waveform inversion.Purely axisymmetric approximation is difficult in practice, however, because the structure along the measurement line of the seismic survey is rarely symmetric with respect to the source location.
In other words, the approach cannot model seismic wave propagation on both sides of the symmetric axis through the seismic source on the measurement line.Furthermore, when one assigns lateral heterogeneity on one side of the cross section, a structural ghost appears on the opposite side because of axisymmetry, such that synthetic seismograms on the side defined as a computation target are contaminated by artificial waves reflected from the ghost that In seismic exploration, treatment of an arbitrary heterogeneous structure model about the axis through a seismic source is crucial for precise comparison between synthetic and observed seismograms, since lateral heterogeneities are close to the axis, and the waveforms calculated by axisymmetric modeling are easily contaminated by artificial reflections from the structural ghost in such a situation.To overcome this difficulty, Takenaka et al. (2003a) proposed a "quasi-cylindrical approach" for seismic exploration, using a nearly linear survey with measurement lines including the source and receiver.In contrast to the conventional axisymmetric approach in cylindrical coordinates using the usual cylindrical domain (0 ≤ r < ∞, −π ≤ φ ≤ π, −∞ < z < ∞), the quasi-cylindrical approach uses a newly defined "quasi-cylindrical domain" Although both approaches calculate the 3-D seismic wavefield on a cross section with a source and receivers, assuming a structure that is invariant in the transverse (φ) direction, the cross section representations are different.In a conventional cylindrical domain, we first have a rectangular half plane with infinite sides formed by movement inside an area specified by ranges 0 ≤ r < ∞ and −∞ < z < ∞, then rotation of this plane in the φ direction through 2π for coverage of the entire spatial domain.Thus, a cross section along the linear survey line of a 3-D target structure is described by two rectangular half planes located at φ = 0 and φ = π.When we assign a 2-D structure model on the φ = 0 plane, the structure on the φ = π plane becomes symmetric, because of the calculation based on axisymmetric modeling.In this situation, the r direction becomes opposite when crossing over the symmetry axis, which makes it impossible to calculate r derivatives in the elastodynamic equation, and therefore the waves cannot travel through the symmetry axis.In fact, conventional axisymmetric modeling produces artificial reflection at the axis, because the line acts rigidly.Nevertheless, such reflection can be regarded as waves coming from the opposite side through the axis in so far as we treat them as axisymmetric wavefields.This is the reason why conventional axisymmetric modeling in cylindrical coordinates cannot treat asymmetric structures with respect to the source axis.On the other hand, the quasi-cylindrical domain first has a rectangular plane with infinite sides formed by −∞ < r < ∞ and −∞ < z < ∞, and then rotates this plane in the φ direction through π to cover the whole domain.In this domain, a cross section of the structure model along the survey line is described by only one plane for φ = 0, and the direction of the horizontal coordinate (r) is unchanged across the vertical axis r = 0 (Figure 3).Hence, we can assign an arbitrary structural model on this plane, followed by reproduction of seismic wavefields propagating through the axis, calculating the r derivatives.The quasi-cylindrical approach, therefore, can calculate realistic 3-D seismic wavefields for an arbitrary cross section of a 3-D structural model with lateral heterogeneity, maintaining the efficiency of conventional axisymmetric modeling.If the structure is defined in a 2-D Cartesian domain (x, z) with shot position x = x 0 , equations for cylindrical coordinates are solved by setting r = x − x 0 .
When we need synthetic seismograms for another shot position in the same structure, we only shift the source grid position in the numerical code, without remaking the computational structure model.Takenaka et al. (2003a) applied the method to a realistic structure model of the Nankai trough off Japan, producing possible observed seismograms by onshore-offshore seismic experimentation in the area.Toyokuni et al. (2005) applied the quasi-cylindrical approach to spherical coordinates.The elastodynamic equation in spherical coordinates is usually solved in the conventional spherical domain (0 ≤ r < ∞,0≤ θ ≤ π, −π ≤ φ ≤ π).However, they introduced a new domain, designated a "quasi-spherical domain" (0 ≤ r < ∞, −π ≤ θ ≤ π, −π/2 ≤ φ ≤ π/2), which maps the sphere in an alternate way.In a conventional spherical domain, we first have a semi-circle with infinite radius formed by rotation from θ = 0t oθ = π, then rotation of this semi-circle in the φ direction through 2π to cover the entire spatial domain.Thus, a cross section along a great circle of the Earth is described by two semi-circles located at φ = 0 and φ = π.When we assign a 2-D structure model on the φ = 0 plane, the structure on a φ = π plane becomes symmetric because of axisymmetry.Similar to the cylindrical case, we cannot take θ derivatives on the source axis, which makes it impossible to propagate waves across the line.On the other hand, in the quasi-spherical domain, we first have a circle with an infinite radius formed by rotation from θ = −π to θ = π, and then rotate this circle in the φ direction through π to cover the entire domain.In this new domain, a cross section along a great circle of the Earth is described by only one circle for φ = 0, and the θ direction is unchanged across the 96 Seismic Waves, Research and Analysis www.intechopen.comsource axis θ = 0 (Figure 4).Hence, we can apply an arbitrary structure model on this plane.We further explain the concept of the quasi-spherical approach in an intuitive way. Figure 4 appears as though the quasi-spherical domain is made by gluing two hemispheres together, although such a joint does not really exist.The reality of the quasi-spherical approach is that it calculates seismic wavefields in two axisymmetric spherical structures, connecting the wavefields only on the axes with θ = 0 and θ = π, which makes us approximately treat wave propagation through an asymmetric structure.In Figure 5, we define a blue semi-circle as structure A and a red semi-circle as structure B. Wavefields propagating in structures A and B are solutions of the elastodynamic equation for axisymmetric structures A and B, respectively.However, computation of θ derivatives at the source axis connects and exchanges wavefields for both structures, which results in apparent treatment of realistic wavefields propagating in an arbitrary asymmetric structure made by combining two semi-circles.This concept is easy to understand with reference to the Riemann surface.When we consider a double-valued function, a two-sheeted Riemann surface should be defined.This surface is made by joining the two sheets crosswise along the "branch cuts", so that values can move from the upper to lower images.Although each sheet is continuous through the branch cuts and the function could have values even on the lower sheet, from the top view of the surface it appears like two sheets glued together.

Spherical coordinates
We call the method of solving the elastodynamic equation in spherical coordinates in the quasi-spherical domain the "quasi-spherical approach".This approach enables modeling of seismic wave propagation in a 2-D slice of a global Earth model of arbitrary lateral heterogeneity, with similar computation time and storage as 2-D modeling, but with full consideration of 3-D wave propagation.Using a method to implement arbitrary moment tensor point sources for conventional axisymmetric modeling (Toyokuni & Takenaka,Fig. 5. Schematic drawings of the concept of quasi-spherical approach.(Left) The elastodynamic equation in spherical coordinates is solved separately for both blue and red axisymmetric structures.However, when both wavefields are connected at the source axis, the resulting wavefield appears to propagate through an asymmetric Earth model.(Right) The concept is similar to the Riemann surface.Viewing from the top, the structure looks as if it is made up by gluing two structures, although these structures are continuous across the joint.

97
Quasi-Axisymmetric Finite-Difference Method for Realistic Modeling of Regional and Global Seismic Wavefield -Review and Applicationwww.intechopen.com2006a), Toyokuni & Takenaka (2006b) simulated the seismic wavefield from the 1994 deep Fiji earthquake.This was done to investigate waveform characteristics observed at Antarctica, propagated through an asymmetric structure with anomalous density and seismic wavespeeds below New Zealand.Toyokuni & Takenaka (2011) extended the quasi-spherical FDM scheme to treat attenuative structures and the Earth's center.Although quasi-axisymmetric modeling can be applied to variety of numerical methods, all previous works developed numerical schemes based on the FDM with second-order accuracy in time and fourth-order accuracy in space, with a staggered-grid formulation.Takenaka et al. (2003a) constructed a staggered-grid scheme for rectangular grids of uniform spacing, for quasi-cylindrical computations of P-SV waves from an explosive source.They used a grid configuration with grid points for v z and normal stress components σ rr , σ φφ , and σ zz located on the axis r = 0, as shown in Figure 6.On the other hand, Toyokuni et al. (2005) and Toyokuni & Takenaka (2006b) used a staggered-grid scheme in spherical coordinates for quasi-spherical computations using nonuniform (Pitarka, 1999) and uniform grid configurations for the vertical (r) and the angular (θ) directions, respectively.Such grid configurations were chosen with smaller vertical grid spacings near interfaces with high contrast of material parameters, e.g., the free surface and the CMB.However, the structural models in these computations were defined over an area with maximum depth 5321 km, so the computations did not treat waves propagating through the Earth center because of problems in this region.The FDM computations of seismic wavefields in spherical coordinates with uniform gridding in the θ direction fail near the Earth center, because of two reasons: (1)

FDM implementation
The extremely small lateral grid spacings near the center perturb the FDM stability criterion, and (2) the singularity of the elastodynamic equation at the center r = 0. To solve the first problem, Toyokuni & Takenaka (2011) applied the so-called multi-domain technique (e.g., Aoi & Fujiwara, 1999;Thomas et al., 2000;Wang & Takenaka, 2001;2010), in which several domains consisting of FD grids with different lateral grid spacings are connected in the r 98 Seismic Waves, Research and Analysis www.intechopen.comdirection, with coarser lateral grids around the center.The second problem was solved using linear interpolation of wavefield variables in the r direction, giving values of particle velocity and stress at the center.Further, Toyokuni & Takenaka (2011) introduced anelastic attenuation into the quasi-spherical FDM.The anelastic behavior of Earth material can be approximated by viscoelastic models, in which the stress-strain relations contain the convolution integral in the time domain, so that time-domain computation such as the FDM had difficulty treating the integral.However, a method using so-called memory variables, which replace the convolution integral with ordinary differential equations for additional internal variables, was invented in 1980s following improvements in Cartesian coordinates (e.g., Carcione et al., 1988a;b;Emmerich & Korn, 1987;JafarGandomi & Takenaka, 2007).Toyokuni & Takenaka (2011) applied the scheme for the first time to the FDM computations in spherical coordinates.The studies with the quasi-spherical FDM used a grid configuration with grid points for v r , σ rr , σ θθ , σ φφ , and σ θφ located on the axis θ = 0, as shown in Figure 7.As mentioned in the previous section, in the staggered-grid scheme, the derivatives of a field quantity are naturally defined halfway between the grid points where the field quantity is defined.Thus, terms on the right-hand side of the elastodynamic equation, including spatial derivatives, are consistently evaluated at the same grid position where the field quantity on the left-hand side is defined.However, this is not the case for terms that do not include spatial derivatives, so these terms have sometimes been evaluated using linear interpolation, despite a decline in accuracy of these terms to second order.To retain fourth-order computation in space at nearly all grid points except along and near the source axis and several computational boundaries, the quasi-axisymmetric schemes prepare the elastodynamic equation that has been rewritten through identities, such as Quasi-Axisymmetric Finite-Difference Method for Realistic Modeling of Regional and Global Seismic Wavefield -Review and Applicationfollowed by discretization (e.g., Takenaka et al., 2003a;Toyokuni et al., 2005).Finally and for example, equations for the P-SV waves corresponding to Eqs. ( 17)-( 22) used by the quasi-cylindrical computations become Similarly, for the quasi-spherical approach, equations for P-SV waves corresponding to Eqs. ( 33)-( 38) can be rewritten as As mentioned above, a characteristic of quasi-axisymmetric modeling is the computation of seismic wavefields even on the source axis.This permits waves to propagate across the axis, from a structure assigned on the right half, to that on the left half of the cross section, and vice versa.For direct computation of the elastodynamic equation on the source axis, we must also solve singularity problems associated with the axis.The elastodynamic equation in cylindrical coordinates has terms containing σ θθ /r and v r /r, which cannot be directly calculated on the axis r = 0. Takenaka et al. (2003a) exploited the formulae derived from limiting operations using the l'Hospital rule, which is also used in Toyokuni et al. (2005) and associated works.For example, formulae for evaluation of wavefield variables on the θ = 0 and θ = ±π axes in quasi-spherical computations are where the variable a can be replaced by σ rθ , σ θθ , σ φφ , σ θφ , v θ ,orv φ .
Since the FDM calculates seismic wavefields only on grid points distributed across computation space, accurate treatment of material discontinuities inside the grid cells has been a serious problem.One possible solution to this problem is the introduction of so-called effective parameters for the density and elastic moduli, calculated by volume arithmetic averaging of densities and volume harmonic averaging of elastic moduli in the cells.The effective parameters scheme enables us to place a material discontinuity at an arbitrary position inside a grid cell (e.g, Boore, 1972;Moczo et al., 2002).Toyokuni & Takenaka (2009) extended the scheme to spherical coordinates and developed a FORTRAN subroutine ACE that calculates the effective parameters analytically for an arbitrary spatial grid distribution within the four major, standard Earth models.Quasi-Axisymmetric Finite-Difference Method for Realistic Modeling of Regional and Global Seismic Wavefield -Review and Applicationwith subducting slab (Takenaka et al., 2003b;c).Figure 8 indicates the P-wave velocity model for the Nankai trough region, Japan, where the Philippine Sea Plate is subducting toward the Eurasian Plate (Kodaira et al., 2002).Each layer in the model has a constant P-wave velocity, corresponding to the color scale.The V P /V S was assumed to be 1.73, except sea water where V P and V S were set to be 1.5 km/s and zero, respectively.Densities for the solid layers were evaluated using the formula of Darbyshire at al. (2000).The model was defined on a 7100 × 1000 grid of spacing 50 m in both horizontal and vertical directions.The time increment was 2.5 × 10 −3 s.We calculated synthetic seismograms for three horizontal positions (S 1 , S 2 , and S 3 ) of 100 m deep seismic sources, as shown in Figure 8.Note that S 1 and S 3 are land and sea shots, respectively, and S 2 is located at the land-sea boundary.The source time function was a phaseless, bell-shaped pulse of width 0.5 s. Figure 9 shows synthetic seismograms for both vertical and horizontal components of particle velocity on the land surface and sea bottom for the three shots.The FDM computations simulated all possible seismic phases in the computation time window.Because of the completeness of the FDM seismograms, we can perform a direct comparison with observed seismograms, which is very important for testing and improving the structural models obtained by seismic surveys.Next, we apply the quasi-spherical FDM for the spherically symmetric Earth model PREM (Dziewonski & Anderson, 1981), to show SH wave propagation for a shear source.Simulation of SH waves is useful for extracting effects related to S waves, since such sources do not exist in nature.It also enables computations for higher frequencies, through a reduction of computational requirements compared to P-SV or 3-D wave simulations.These attributes are why many authors have been working with global SH-wave computations (e.g., Chaljub & 104 Seismic Waves, Research and Analysis Tarantola, 1997;Igel & Weber, 1995;Igel & Gudmundsson, 1997;Jahnke et al., 2008;Thorne et al., 2007;Wang & Takenaka, 2011;Wysession & Shore, 1994).We use a shear point source M 21 = −M 12 at depth 600 km, with source time function described as a phaseless, bell-shaped pulse of width 4 s.The computational model is defined on a 2843(r) × 27000(θ) grid with maximum depth 2891 km (CMB), since SH waves cannot propagate into the outer core.The free surface condition has been applied to the top and bottom boundaries of the computational domain.Uniform grid spacing is used in both vertical and angular directions.The time increment is 7.0 × 10 −3 s. Figure 10 shows sequential snapshots at six time steps, which allows us to confirm fundamental properties of SH-waves reverberating inside the crust and mantle.

Applications
Finally, to investigate characteristics of observable waveforms in the intra-Antarctic region, we calculate global synthetic seismograms with the quasi-spherical FDM for an asymmetric model with a simply-shaped, high seismic wavespeed anomaly superimposed on the attenuative PREM.Numerous temporal broadband seismic stations have been recently installed in this region, in association with the International Polar Year (IPY) 2007-2008 (Kanao et al., 2009).The seismic source is located in northern Bolivia at depth 651 km, the same location as the 1994 deep Bolivia earthquake.However, the mechanism is a simple dip-slip source with nonzero moment-tensor components of M 13 = M 31 .The source time function is a phaseless, bell-shaped pulse with duration 60 s.The anomaly is expressed as a region containing perturbations on P-and S-wavespeeds set at +20.0 % above the PREM basis, within vertical and angular ranges of 3480km ≤ r ≤ 4180km and 16.18 • ≤ θ ≤ 46.18 • .This is representative of a high velocity anomaly beneath southern South America, deduced by seismic tomography.We calculate wavefields on a longitudinal cross section, including the source and the anomaly, to see how the observed seismograms in Antarctica are affected by the anomaly.As shown in Figure 11, the angular range of Antarctica for this situation is 52.78 • ≤ θ ≤ 99.58 • .We use the ACE subroutine (Toyokuni & Takenaka, 2009) to generate the effective parameters for the PREM, with respect to the given grid distribution in the radial direction.The results are presented by synthetic seismograms along the Earth surface, and sequential snapshots of wave propagation.Figures 12 and 13 indicate respectively the angular (θ) and transverse (φ) components of synthetics at various angular ranges in Antarctica, calculated for (a) the PREM and (b) the model, with a high velocity anomaly.Differential seismograms in panel (c) are obtained by subtracting the PREM results from those of the asymmetric model, which clearly illustrate various phases affected by the anomaly region.Since the anomaly is located just above the CMB, we see that the core reflection such as ScS, sScS, and their multiple reflections, have been strongly affected by the region, as expected.These results suggest probable characteristics of observed seismograms in the intra-Antarctic region.Figure 14 shows sequential snapshots of the vertical (r) component of the seismic wavefield propagating on a cross section at every 300 s, from 300 s to 3900 s after excitation.We can see the asymmetric wavefield about the source axis, caused by the anomaly.The computation required 2.4 Gbytes of memory in a single precision calculation, with computation time of 27.3 hours on eight CPUs with IBM POWER6 architecture (4.7 GHz clock speed), for a total duration of 5000 s after excitation.

105
Quasi-Axisymmetric Finite-Difference Method for Realistic Modeling of Regional and Global Seismic Wavefield -Review and Application -

Conclusions
We have reviewed recent developments of numerical computation methods that have accuracy and computational efficiency for realistic seismic wavefields, using the FDM.
Traditional axisymmetric modeling solves the 3-D seismic wave propagation only on a 2-D cross section of a structure model including a seismic source and receivers, under the assumption that the structure is invariant in the transverse direction about the axis through the source.However, realistic structures with asymmetry cannot be treated in principle.
Quasi-axisymmetric modeling represents methods solving the seismic wave equation in newly defined quasi-cylindrical / spherical coordinates, rather than the usual cylindrical / spherical coordinates.This type of modeling retains the efficiency of axisymmetric modeling but can treat an arbitrary asymmetric structure, thereby providing a breakthrough for the problem of traditional axisymmetric strategies.

Fig. 6 .
Fig.6.Staggered-grid distribution used in quasi-cylindrical FDM computations ofTakenaka et al. (2003a).The grids for the vertical component of particle velocity v z and the normal stress components σ rr , σ φφ , σ zz are located on the source axis r = 0. ∆r and ∆z are grid spacings in the r and z directions, respectively.

Fig. 7 .
Fig.7.Staggered-grid distribution used in quasi-spherical FDM computations of, for example,Toyokuni et al. (2005).The grids for the vertical component of particle velocity v r , the normal stress components σ rr , σ θθ , σ φφ , and the {θφ}-component of the stress tensor σ θφ are located on the source axis θ = 0. ∆r and ∆θ are grid spacings in the r and θ directions, respectively.

Fig. 8 .
Fig. 8. P-wave velocity model used for simulation of onshore-offshore seismic experiment.Stars indicate shot locations (after Takenaka et al., 2003b).This section shows examples of wavefield computation using quasi-axisymmetric modeling.First, we display an application of the quasi-cylindrical FDM to a realistic structure model

Fig. 10 .
Fig. 10.SH-wave snapshots up to a period of4satsixtime steps, showing generation and propagation of various seismic phases in the spherically symmetric Earth model PREM.Each frame uses the same color scale: red and blue indicate plus and minus amplitudes, respectively.Solid circles are the free surface, the 670-km discontinuity, and the CMB.Seismic source is a 600-km deep shear point source, indicated by a star. 103

Fig. 11 .
Fig. 11.Cross section of structural model for computation of synthetic seismograms.Circles indicate the free surface, 400-km and 670-km discontinuities, the CMB and ICB.Seismic source is located at depth 651 km, underneath northern Bolivia.Blue area is an anomaly of seismic wave speed, placed just above the CMB within a range of 3480km ≤ r ≤ 4180km and 16.18 • ≤ θ ≤ 46.18 • , having a +20% velocity increase of P-and S-wavespeeds from the PREM basis.A red line indicates angular range of Antarctica (52.78 • ≤ θ ≤ 99.58 • ).

Fig. 12 .
Fig. 12. Synthetic seismograms of v θ at the Earth surface within the angular range of Antarctica, calculated for (a) the PREM, and (b) the model, including a high velocity anomaly.Differential seismograms (c) are calculated by subtracting (a) from (b), which indicate various phases affected by the anomaly.All traces were low-pass filtered with cutoff period 60 s, and are shown at the same scale.

Fig. 13 .
Fig. 13.Synthetic seismograms of v φ at the Earth surface within the angular range of Antarctica, calculated for (a) the PREM, and (b) the model, including a high velocity anomaly.Differential seismograms (c) are calculated by subtracting (a) from (b), which indicate various phases affected by the anomaly.All traces were low-pass filtered with cutoff period 60 s, and are shown in the same scale.