Modelling seismic wave propagation in a two dimensional cylindrical whole earth model using the pseudospectral method

SUMMARY

We present a method for modelling seismic wave propagation in a whole earth model by solving the elastodynamic equations in 2 D cylindrical coordinates (r, θ) using the Fourier pseudospectral method (PSM). In solving the 2 D cylindrical elastodynamic equations for a whole earth model, a singularity arises at the centre (r=0) of the earth. To avoid the singularity, we develop a scheme that uses extension of field variables in the radial direction, with which computation of the wavefield at the centre is avoided, so that the wave propagation through the centre can be calculated. The time interval used in the calculation is determined by the smallest lateral grid spacing around the centre in the model. In a cylindrical coordinate system, the smallest lateral grid spacing is generally so small that the calculation is too time consuming to be realistically carried out even on a supercomputer. We adopt a multidomain scheme to increase the smallest lateral grid spacing and avoid the oversampling of the physical domain around the centre of the earth. A smoothing scheme in the wavenumber domain is also proposed, which enables us to use a large enough time interval to allow the calculation for the whole earth model on a desktop workstation. The waveforms calculated by the present method are compared with those obtained by the Direct Solution Method (DSM) to demonstrate their high accuracy. This method significantly reduces the computer memory and computation time required and makes it possible to study the effects of small wavelength heterogeneities that can be approximated as azimuthally symmetric on wave propagation in the earth. We apply the present method to study the effects of local heterogeneity in the earth by adding a low velocity perturbation above the core–mantle boundary (CMB) to the IASP91 earth model.

1 Introduction

With the accumulation of high quality broad band global seismic data, recent studies of whole earth structure have revealed heterogeneity in both the radial and lateral directions in the mantle (e.g. Hedlin et al. 1997) and at the core–mantle boundary (CMB) (e.g. Lay et al. 1998), and anisotropy in the inner core (e.g. Song 1997). Such discoveries have been made from the detailed analysis of observed waveforms for a variety of phases that propagate through, are reflected, refracted, converted or diffracted due to the heterogeneity, and the associated anomalies in traveltimes and amplitudes. Forward modelling of seismic wave propagation in the whole earth with arbitrary heterogeneity and the synthetic seismograms at the ground surface are crucial in the verification and understanding of the observed phase anomalies, further constraining the heterogeneity derived from observations.

Modelling seismic wave propagation in the spherical whole earth has been carried out using several methods (Alterman et al. 1970; Li & Tanimoto 1993; Wysession & Shore 1994; Friederich & Dalkolmo 1995). With developments in computer technology and numerical simulation algorithms, recent modelling has been performed for spherical whole earth models with arbitrary heterogeneity of strong velocity perturbations. For instance, Yoon & McMechan (1995) applied the finite difference algorithm on a staggered grid in Cartesian coordinates to simulate the wave propagation inside a 3 D whole earth model. They calculated the complete wavefield for a PREM model with a bump on the CMB, to illustrate the effects of the CMB topography on synthetic seismograms for long period (80 s) body waves. However, the accuracy of their calculations was not checked by comparison with those obtained by other methods. Igel & Weber (1995, 1996) performed simulations of seismic wave propagation in the mantle for axisymmetric earth models using a finite difference scheme on a staggered grid in spherical coordinates. They studied the effect of CMB structure on the seismograms by comparing three different velocity models with lowermost mantle heterogeneity. Igel & Gudmundsson (1997) extended their method to simulate SH wave propagation in the mantle with a depth dependent lateral grid spacing. They studied the frequency dependence of arrival times of long period S and SS waves through random upper mantle models with specified spectral properties. Chaljub & Tarantola (1997) used the staggered grid finite difference scheme to study the topography effect of the upper mantle 660 km discontinuity on SS precursors. However, the schemes presented by Igel & Weber (1995, 1996) and Chaljub & Tarantola (1997) were only applied to wave propagation in the mantle, and have not been used to simulate wave propagation through the core of the earth. Since the lateral grid spacing decreases with depth for models defined in spherical coordinates, there is a maximum limit to the depth range that can be simulated for a reasonable time interval. Cummins et al. (1997) applied the Direct Solution Method (DSM) (Geller & Ohminato 1994; Cummins et al. 1994a,b; Takeuchi et al. 1996) to a laterally heterogeneous spherical earth model with strong axisymmetric velocity perturbations in the upper mantle of the IASP91 model, and showed the effect of the heterogeneity on S wave propagation in the mantle. The Spectral Element Method (SEM) (e.g. Komatitsch & Vilotte 1998; Komatitsch & Tromp 1999), which is very flexible in handling free surface topography, fluid–solid interfaces, anisotropy and attenuation with high accuracy has been applied to wave propagation simulations in 3 D whole earth models (e.g. Chaljub & Vilotte 1998; Capdeville et al. 1999). The Chebyshev spectral method was used by Igel (1999) to simulate seismic wave propagation in the 3 D spherical earth with heterogeneity in the uppermost mantle defined in spherical coordinates, but the range of the physical domain was limited to 80° in both the radial and lateral coordinates and 5000 km in depth. Furumura et al. (1998) exploited a pseudospectral method (PSM) scheme to simulate seismic wave propagation in laterally heterogeneous whole earth models. They solved the elastodynamic equations for a 2 D heterogeneous earth structure in 2 D cylindrical coordinates whose physical domain ranged from the earth's surface to 5315 km below the surface, including the outer core. They applied the method to the IASP91 model and predicted nearly all of the seismic phases in the whole wavefield along the surface of the earth. They also applied their method to 2 D heterogeneous models to study the anomalies in arrival times and amplitudes of various phases caused by heterogeneity located just below the earth's surface and above the CMB. Since a 2 D approximation was made, they could simulate the whole wavefield for relatively short period (15 s) in the entire cross section through a great circle of the earth except the core. This 2 D cylindrical method was used by Furumura et al. (1999) to study the effect of stochastic heterogeneity on seismic wave behaviour in the mantle. They compared the effects of broad scale and stochastic heterogeneity for a model built from a slice through a tomographic model for the Himalayan region. However, in this 2 D cylindrical method, the inner core cannot be included in the model because of the singularity at the centre of the earth and the limitation in the minimum grid spacing related to the time interval. Therefore, the seismic wave propagating through the inner core could not be calculated. Thomas et al. (2000) recently presented a scheme to solve the ‘acoustic’ wave equation in spherical coordinates for axisymmetric media using a high order finite difference method. They developed a multidomain approach to avoid the stability problem associated with the centre of the earth, and proposed a scheme to treat the centre of the earth in a Cartesian system that allows the wave propagation through the whole core to be calculated. This method was applied to an axisymmetric earth model to simulate P wave propagation in the whole earth and the effects of scatterers in the earth's lower mantle on core phases.

In this paper, we present a method to simulate seismic wave propagation in a heterogeneous whole earth model including the core using the Fourier pseudospectral method (Kosloff & Baysal 1982; Kosloff et al. 1984). We define the model and solve the elastodynamic equations in 2 D cylindrical coordinates. The model is defined between 0 and 2π in the lateral (θ) direction and between the centre of the earth and the surface in the radial (r) direction. We compare the synthetic seismograms calculated by our method with those obtained with the DSM method. If the results of our method are proved to be close enough to the results given by the DSM method, we think it would be reasonable to perform the modelling of the effect of random heterogeneity in the whole earth using a 2 D cylindrical model instead of a 3 D spherical model. Since use of a 2 D model will significantly reduce the required computer memory and computation time, simulations could be performed for models with finer grid spacing for relatively shorter wavelength seismic waves than the full 3 D calculation. However, such a model cannot accurately predict the scattering and focusing/defocusing effects from out of plane wavefields.

To treat the whole earth model including the core, we meet two challenges: one of them is the singularity at the centre of the earth where r equals zero; the other is the reduced lateral grid spacing close to the centre, which leads to a very small time interval for stable computation, which will cause very long computation times for synthetic seismograms of reasonable duration. For the first challenge, we develop an ‘extension scheme’ along the radial direction in which the spatial variables are extended and the spatial derivatives are calculated along the range of the diameter instead of the radius. For the second one, we adopt a multidomain scheme in which the lateral grid spacing varies with depth and a smoothing technique in the wavenumber domain is also applied for a small region around the centre when performing the spatial derivatives in the lateral direction.

The spatial derivatives in both the radial and lateral directions are calculated in the wavenumber domain by the Fast Fourier Transform (FFT). Compared with the traditional finite difference scheme, the Fourier differential operator can achieve results of the same accuracy with larger grid spacing, so that the computational memory and time can be considerably reduced (e.g. Fornberg 1987; Daudt et al. 1989; Vidale 1990). Kessler & Kosloff (1990, 1991) presented the pseudospectral method for solving the acoustic wave equation and the elastodynamic equation in 2 D cylindrical coordinates. In their method, the derivatives along the lateral direction are calculated by the Fourier expansion, while the derivatives along the radial direction are calculated by the Chebyshev expansion because it can represent the free surface accurately. Since the design of a Chebyshev mesh that is suitable for both computational requirements and the representation of the structural heterogeneity in a whole earth model is rather difficult (Furumura et al. 1998), we use the Fourier expansion for the radial derivatives. The accurate representation of the earth's free surface and the discontinuities in the whole earth model are accomplished by a mapping technique (Fornberg 1988; Furumura et al. 1998) along the radial direction.

In the following sections, we will first explain the scheme for the whole earth modelling, the treatment of the singularity at the earth's centre, the multidomain approach and the smoothing scheme. Next we check the accuracy of the method by comparison of the synthetic seismograms with those obtained by the DSM. We then show some examples in which we will apply the method first for the IASP91 earth model to give a complete image of seismic wave propagation in the whole earth and then for models with strong heterogeneity located above the CMB to see its effect on the whole wavefield.

2 Equations of motion

We consider P–SV wave propagation in a whole earth model. In a cylindrical coordinate system with coordinates (r, θ, z), assuming invariance in z for all fields leads to the 2 D cylindrical equations of momentum conservation for P–SV waves (Aki & Richards 1980),

where ü_p=ü_p(r, θ, t) (p=r, θ) are the acceleration in the radial (r) and lateral (θ) directions at a gridpoint (r, θ) at time t, ρ=ρ(r, θ) is the mass density, f_p=f_p(r, θ, t) (p=r, θ) are body forces and σ_pq=σ_pq(r, θ, t) (p, q=r, θ) are the stress components. The constitutive relations between the stress and the displacement for an isotropic linear elastic solid are

where u_p=u_p(r, θ, t) (p=r, θ) are the displacement components and λ=λ(r, θ) and μ=μ(r, θ) are the Lamé constants.

3 Numerical implementation in 2 D cylindrical coordinates

3.1 Physical and numerical domains

We consider a 2 D heterogeneous earth structure defined between r=0 and 6371 km and θ=0° and 360° that is a slice cutting through the great circle of a 3 D spherical earth. To solve the above equations for such a model, the physical domain is discretized in the radial (r) and lateral (θ) directions in 2 D cylindrical coordinates as depicted in Fig. 1. In the lateral direction, the field quantities g(mΔθ) (m=0, 1,…, M−1) are distributed with uniform angular interval Δθ for given radius r. As shown in Fig. 1, the numerical model is composed of three subdomains with a different number of gridpoints in the lateral direction, and therefore Δθ varies for each subdomain. The multidomain approach used here will be described in the next subsection. In the radial direction, the quantities g(nΔr) (n=0, 1,…, N−1) are distributed along grids with irregular spacing Δr in depth that is small at the free surface and interfaces in the model. The varying radial grid spacing Δr is achieved by using a mapping technique (Fornberg 1988; Furumura et al. 1998) in order to locate accurately the free surface and the interfaces in the earth. The actual radial grid spacing Δr and the lateral grid arc length along the radius are shown together with the model in Fig. 1.

3.2 Multidomain approach

In a cylindrical coordinate system, with decreasing radius from the surface to the centre, as can be seen in Fig. 1, the lateral grid spacing (the arc length between adjacent gridpoints) decreases very quickly since it is calculated by 2πr/M, where M is the number of gridpoints in the lateral direction. In the case of a single numerical domain, in which the number of lateral gridpoints is the same for radius r, the physical area covered by the lateral gridpoints becomes smaller from the surface to the centre but the number of gridpoints remains the same. This leads to oversampling of the physical domain near the centre of the earth, and such oversampling will occupy much unnecessary computer memory and cause extra execution time in the calculation.

The time interval determined by the stability condition is selected according to the smallest grid spacing in the numerical mesh. In cylindrical coordinates, the smallest grid spacing is the lateral one at gridpoints nearest to the centre. In the whole earth model including the core, radius r of the first gridpoint next to the centre should be small compared to the minimum wavelength in the modelling. This will cause extremely small lateral grid spacing next to the centre and the related time interval will then be extremely short, which means the modelling cannot realistically be performed, as will be shown in the following examples. Therefore, increasing the smallest lateral grid spacing is crucial if the earth's core is included in the model.

In order to treat the problems mentioned above, we adopt a multidomain approach. The whole numerical domain is divided into several subdomains with different lateral grid spacing as shown in Fig. 1. The subdomain around the centre has fewer lateral gridpoints and larger grid spacing compared with the subdomain near to the surface. In the examples given in the following sections, the total computational domain is composed of three subdomains. The numbers of lateral gridpoints are 256, 512 and 1024 for subdomains from the centre to the surface. The lateral grid spacing decreases by half between two adjoining subdomains. In this case, the smallest lateral grid spacing will increase fourfold and the time interval will also increase fourfold, and the computation time will decrease fourfold compared with the model in which only a single domain (1024 lateral gridpoints) is used. This multidomain method also samples the physical domain evenly for all the area as shown in Fig. 1. Kessler & Kosloff (1991) used a multidomain technique in simulations of elastic wave propagation in the vicinity of cylindrical objects in 2 D cylindrical coordinates. In their study, the radial derivatives are solved by the Chebyshev expansion and the adjacent subdomains are connected by the characteristic variables of the wave equation (that is, the components of the displacement and the stress). In our method, the discontinuities of the field variables across the boundary between two subdomains are connected through interpolation of the field variables performed using the FFT. We interpolate the field variables in the subdomains with fewer lateral gridpoints to the same number as the grid number of the outermost subdomain and then calculate the radial and lateral derivatives. Since the number of lateral gridpoints in each subdomain is a power of 2 defined between 0 and 2π in the PSM method, the interpolation using FFT is quite easily incorporated in the scheme.

3.3 Avoiding the singularity at the centre

Eqs (1) and (2) include terms divided by r, therefore a singularity arises at the centre where r equals zero when we consider the whole earth model including the core. We exploit the following extension of field variables to avoid this singularity and calculate the radial derivatives.

We observe the field variables defined along one diameter in the model at θ and θ+π as shown in Fig. 2. Let g(r, θ) represent the field variable defined along the radial direction in cylindrical coordinates, where g(r, θ) is defined for 0≤r≤R and 0≤θ≤2π, where R is the range of the physical domain. In order to calculate the wave propagating through the centre, we solve the wave equations in the range −R≤r≤R instead of 0≤r≤R. To do this, we map g(r, θ) defined along the radius for 0≤r≤R and 0≤θ≤2π to g′(r′, θ′) defined along the diameter for −R≤r′≤R and 0≤θ′≤π. We write the discrete form of g(r, θ) along the radius as g_θ(i), and consider g_θ(i) and g_θ+π(i) (i=0, 1,…, N−1) distributed on the diameter (Fig. 2). For the mapping, we reverse the order of g_θ+π(i) for 0≤r≤R at θ+π to form the first half (corresponding to −R≤r′<0) of a new array g′(j)(j=0, 1,…, 2N−1) at θ′=θ, while g(i) for 0≤r≤R at θ are put directly into the second half (0<r′≤R) of g′(j), as shown in Fig. 2.

The radial derivatives are then evaluated for the new array g′(j) in the diameter range (−R≤r′≤R, 0≤θ′≤π). The signs of the radial derivative operation and some field variables are changed on the mapping for the range 0≤r≤R at θ+π as listed in Table 1. The signs of u_r and σ_rθ are changed because the positive direction of r is reversed, while the signs of u_θ and σ_rr remain unchanged. The sign of the radial derivative operation ∂/∂r is changed due to the reversed order of the field variables.

After calculating the radial derivatives, the mapping is again performed for the ∂/∂r terms on −R≤r′≤0 at θ′=θ because the following calculations using these derivatives are carried out for the original range 0≤r≤R and 0≤θ≤2π. The changes of the signs of the radial derivatives for −R≤r′≤0 at θ′=θ are given in Table 2. Since the signs of u_r and σ_rθ changed before taking the derivatives as shown in Table 1, the signs of ∂u_r/∂r and ∂σ_rθ/∂r will not change because of the sign reversal of ∂/∂r. The signs of u_θ and σ_rr are not changed before taking the derivatives so the ∂u_θ/∂r and ∂σ_rr/∂r will change their signs.

3.4 Spatial derivatives and time extrapolation

Both of the spatial derivatives in the radial and lateral directions are calculated by multiplication in the wavenumber domain, and the transformation between the physical domain and the wavenumber domain is performed by the FFT (Kosloff et al. 1984; Furumura et al. 1998). Since we use the multidomain approach and the ‘extension scheme’ to avoid the singularity at the centre, the calculation of the spatial derivatives and the stress and displacement components in eqs (1) and (2) at each time step is accomplished by the following procedures.

• (i)Compute ∂/∂θ of u_r and u_θ in each subdomain for 0≤r≤R, 0≤θ≤2π.

• (ii)Perform interpolation in the θ direction in the subdomains by FFT as mentioned in Section 3.2 to form the array g(i, j)(i=0, 1,…, N−1); (j=0, 1,…, M−1) for u_r and u_θ, where N and M are the number of radial and lateral gridpoints, respectively.

• (iii)Use the ‘extension scheme’ as mentioned in Section 3.3 to map g(i, j) to g′(i′, j′) (i′=0, 1,…, 2N−1); (j′=0, 1,…, M/2−1), where 2N is the number of radial gridpoints for −R≤r′≤R, and M/2 is the number of lateral gridpoints for 0≤θ′≤π.

• (iv)Compute ∂/∂r′ of the mapped u_r′ and u_θ′ for −R≤r′≤R, 0≤θ′≤π. Map the ∂/∂r′ back to ∂/∂r terms again for 0≤r≤R, 0≤θ≤2π.

• (v)Calculate σ_rr, σ_rθ and σ_θθ from u_r, u_θ and their ∂/∂r and ∂/∂θ terms computed in steps (i) and (iv) over 0≤r≤R, 0≤θ≤2π.

• (vi)Repeat steps (i) to (iv) for σ_rr, σ_rθ and σ_θθ.

• (vii)Compute the acceleration ü_r, ü_θ from σ_rr, σ_rθ, σ_θθ and their ∂/∂r and ∂/∂θ terms in eq. (1) for the current time step.

• 3

  (viii)Integrate the acceleration twice to obtain the displacement u_r, u_θ using the following second order finite difference scheme:

  

  

  3

  where p=r, θ, ü_pⁿ are the accelerations calculated in the current time step in (vii) and u_pⁿ and u_pⁿ⁺¹ are displacements at the current and the next time step, respectively.

In the radial direction, the field variables are distributed on irregularly spaced grids. The ∂/∂r terms, which are first calculated for evenly spaced gridpoints by FFT, are converted to the values on irregular spaced grids by the mapping technique used by Fornberg (1988) and Furumura et al. (1998). The ∂/∂θ in the θ direction are calculated in each subdomain with a different number of gridpoints with a uniform grid spacing.

3.5 Seismic source

We calculate the wavefield excited from a point source in 2 D cylindrical coordinates that is equivalent to a line source extending infinitely in the z direction in (r, θ, z) coordinates. The body force corresponding to a combination of moment tensor components M_rr(t), M_rθ(t) and M_θθ(t) is introduced over a small region of the gridpoints around the centre of the source (r₀, θ₀). The body force system for a moment tensor point source in 2 D cylindrical coordinates is derived by differentiating the equations for single point forces in 2 D cylindrical coordinates with respect to the source coordinates (e.g. Ben Menahem & Singh 1981) as follows:

The 1/r and 1/rr₀ can also be replaced by 1/r₀ and 1/r₀² in a distribution sense, because of the following property:

The source expression (4) is different from that used by Furumura et al. (1998, 1999), which did not have the factors 1/r and 1/rr₀. In this paper, the delta functions δ(r) and δ(θ) in eq. (4) are approximated by Herrmann's pseudo delta function (Herrmann 1979), which gives a point source of unit area in (r, θ) coordinates. We also use Herrmann's function for the source time history.

3.6 Boundary conditions

In the whole earth model including the core, the free surface boundary of the entire numerical domain only needs to be considered in the radial direction. The free surface is incorporated into the model by satisfying the zero traction condition (σ_rr=σ_rθ=0) at the surface. We adopt an ‘image method’ to implement this condition, which was proposed by Crase (1990) and successfully employed by Rodrigues & Mora (1993). We add a number of gridpoints above the free surface. σ_rr and σ_rθ at these gridpoints are obtained by using the anti symmetric extension from the values on corresponding gridpoints below the surface, while the displacement components u_r and u_θ are symmetrically extended from the values at corresponding gridpoints to ensure stability (e.g. Crase 1990; Rodrigues & Mora 1993). This scheme is different from the ‘symmetric differentiation’ that Furumura et al. (1998, 1999) employed for the free surface condition. Therefore, the simulation can be performed stably for a long duration of seismograms and the stability problem related to the free surface condition as mentioned by Furumura et al. (1999) does not appear. Robertsson (1996) showed that this image method is very accurate. In this paper, the total number of gridpoints in the radial direction is 256, and 64 of them are located above the free surface. The liquid outer core is treated by setting μ=0. This treatment for a liquid layer has been successfully applied to seismoacoustic scattering problems (e.g. Bayliss et al. 1986; Dougherty & Stephen 1988, 1991; Levander 1988; Stephen & Swift 1994; Swift & Stephen 1994; Okamoto 1994; Robertsson & Levander 1995; Robertsson et al. 1996; Okamoto & Takenaka 1999). We reduce the radial grid spacing at the CMB and inner core boundary (ICB) by mapping, as mentioned in Section 3.1, which is quite efficient in suppressing the Gibbs' noise produced during the differentiation calculation using the FFT at the discontinuities (CMB and ICB), since the oscillation noise decays exponentially as the number of gridpoints from the boundary increases (Furumura et al. 1998). In the θ direction, the periodic boundary condition is naturally incorporated in the calculations from the periodicity in the FFT.

4 Comparison with the Direct Solution Method

In this section, we check the validity and accuracy of our method by comparing our results with those obtained using the DSM (Geller & Ohminato 1994; Cummins et al. 1994b; Takeuchi et al. 1996), which gives exact waveforms for spherically symmetric media. The numerical model we used for the comparison calculation is composed of three subdomains, as shown in Fig. 1. The first subdomain around the centre of the earth extends to the middle of the inner core, and the lateral grid spacing ranges from 0.43 to 13.54 km. The second subdomain includes part of the inner core and most of the outer core, and the lateral grid spacing is between 7.11 and 37.30 km. The third subdomain covers the upper part of the outer core and the mantle and extends to the surface, where the lateral grid spacing varies from 18.85 to 39.09 km. The radial grid spacing changes from 28.10 km (in the inner core, at the CMB and at the surface) to 39.10 km (in most of the mantle and the middle of the outer core). The numbers of gridpoints in the lateral direction for the three subdomains from the core to the surface are 256, 512 and 1024, respectively, and 256 in the radial direction. The source is a double couple point source located at a depth of 600 km with moment tensor components M_rr=−1.0 and M_θθ=1.0 (the other components are all zero). The width of the source time function is 50 s. Using the minimum S wave velocity just below the surface (V_S=3.36 km s⁻¹) and the minimum S wavelength (168.0 km) for the source, the number of gridpoints per minimum wavelength for the largest grid spacing is 4.3. The time interval Δt used in the calculation is constrained by the ratio of the minimum grid spacing in the model and the maximum wave velocity as

The minimum grid spacing is 0.43 km near the centre and the maximum wave propagation speed is V_P=13.7 km s⁻¹, so the time interval will be 0.008 s if we use α=0.26 for a 1 per cent tolerance error level (Daudt et al. 1989). This value of Δt is too small for actual calculations.

In order to increase Δt, we apply smoothing in computing the lateral derivatives. Since the grid spacing in the lateral direction around the centre in the first subdomain is still too small compared with the minimum S wavelength even when the multidomain is used, we apply a low pass filter to filter out the high wavenumber component when we calculate the lateral derivatives in the wavenumber domain. For lateral gridpoints in a circle at distance r from the centre, the Nyquist wavenumber along the circle is

where Δa is the lateral arc length between adjacent grids. We take a reference lateral arc length Δa^ref, and define the cut off wavenumber as

When we calculate the lateral derivatives for lateral grids at r, we filter out the wavenumber component higher than K^c and the time interval Δt in eq. (6) is then determined by Δa^ref instead of Δ^min. Since the Δa^ref is much larger than the Δ^min in the model, the Δt actually used in the calculation will be increased by an order of magnitude. The accuracy of the results is not affected by this smoothing if the Δa^ref is selected to be small enough compared with the minimum wavelength in the modelling.

In the comparison modelling, we set Δa^ref to 14.0 km, which leads to 12 gridpoints per minimum S wavelength. The reference lateral grid spacing then allows a time interval up to 0.27 s. The time interval we used in the computation is 0.25 s, which increases by about 31 times over the time interval (0.008 s) without smoothing. In Fig. 1, the dashed line in lateral grid spacing is the actual arc length in the model and the corresponding solid line is the Δa^ref that is used in the computation to determine the time interval Δt. With 14.0 km as the reference arc length, the smoothing in computing lateral derivatives is actually applied to a very small physical domain within about 570 km from the centre in the inner core. The velocity model employed here is the IASP91 model (Kennett & Engdahl 1991) without considering anelastic attenuation.

Since the DSM solves the equations of motion in 3 D spherical coordinates, for comparison we map the ‘line’ source solution obtained by our method in 2 D cylindrical coordinates to an approximate 3 D ‘point’ source solution (e.g. Vidale et al. 1985; Helmberger & Vidale 1988; Pitarka et al. 1994, 1996). The mapping of the seismograms is performed through the following filter:

where * is the convolution operation, u_PSM(t) is the displacement obtained by our 2 D method for a double couple ‘line’ source, and v(t) is the converted waveform that corresponds to the displacement excited by a double couple ‘point’ source. R is the distance between the source and the observation position, i.e. the hypocentral distance. Furumura et al. (1998, 1999) also compared their 2 D synthetic seismograms with those obtained by the DSM. For the comparison, they corrected only for the difference in geometrical spreading between 2 D and 3 D wave propagation by multiplying by a factor R^−0.5 (R is the epicentral distance in their paper, not hypocentral distance). Their comparisons did not show good agreement. This may be mainly because they did not apply the correction filter for the difference in the pulse shape between the ‘line’ and the ‘point’ source solution [(1/π) ] in eq. (9).

4.1 Comparison of the complete seismograms

We calculated the synthetic seismograms at five epicentral distances (Δ=30°, 60°, 90°, 120° and 150°) by our method, then converted them to 3 D seismograms by applying the filter in eq. (9) and compared the resulting seismograms with the 3 D synthetic seismograms calculated by the DSM. The comparison of the seismograms is shown in Fig. 3. A bandpass filter of 40–300 s has been applied to all the seismograms. In the radial component u_r, we see both seismograms are almost identical in waveforms and traveltimes for all major phases. In the lateral component u_θ, most major phases also show very good agreement in both waveforms and traveltimes, but there is a small phase delay for surface multiples such as the SS phase as observed at Δ=60°, 90° and 120°. This is due to the larger grid spacing in the lateral direction near to the surface in our model. Since the number of gridpoints per minimum S wavelength in the θ direction is 4.3 along the surface, which is less than the six gridpoints per minimum S wavelength in the radial direction, the SS phase that propagates along a shallow zone below the surface is affected by the greater lateral grid spacing just below the surface. The phase delay grows stronger with increasing epicentral distance. However, this could be prevented by using finer grids in the model.

In Fig. 3, large discrepancies for some later phases between the two methods at short epicentral distance are also seen; for example, the P′P′ at around 2342 s in u_r and the (ScS)₂ at about 2000 s in u_θ at 30°, and the P′P′ between 2200 and 2500 s in u_r and u_θ at 60°. The amplitudes of these phases are larger than the results by the DSM. These later phases travel very long distances (e.g. about 25 000 km for P′P′ at 30°) and a long time in the earth, so the accumulated numerical error causes large differences in the amplitudes of these phases.

4.2 Comparison of core phases

Since the whole core of the earth has been included in the model, complete core phases can be calculated by this method. In this subsection, we select some core phases at various epicentral distances from the complete seismograms obtained by the method described above and compare them with the DSM results.

The PKP and pPKP phases that propagate through the inner and outer core can be observed at epicentral distances of 116° to 244°. In Fig. 4(a), the complete PKP and pPKP phases are displayed between 110° and 180°. The agreement in arrival time and amplitude between the two results is fairly good for all epicentral distances, which suggests that these core phases can be correctly calculated using this method. In Fig. 4(b), the P_diff and pP_diff phases diffracted at the CMB are displayed between 110° and 130°. Both the arrival times and the amplitudes of these phases from the two methods coincide with each other very well. Fig. 5(a) shows the PKS and pPKS phases between 120° and 155°, which were converted at the CMB when PK and pPK penetrated from the outer core into the mantle. In Fig. 5(b), the comparison of phases related to the direct S phase are shown between 60° and 115°. The triplication of the direct S, SKS and S_diff occur at around 80°, while SKS and S_diff are observed at epicentral distances larger than 80°. The good agreement for these phases can be seen in the lateral component and for the SP phase in the radial component. As can also be seen in Fig. 3, the surface multiples SS show a delay in arrival time in the PSM results compared with the DSM results.

4.3 Comparison of surface waves

Comparisons of the Rayleigh waves excited from the 600 km deep source in the calculations are shown in Fig. 6(a). The calculations have been performed for seismograms of 3000 s duration, where the surface waves can be seen up to an epicentral distance of 90°. In Fig. 6(a), the Rayleigh waves between 10° and 90° are displayed by superimposing the PSM results onto the DSM results. We see that both the amplitudes and the arrival times of the Rayleigh waves are nearly identical for the two results up to 80°, and show discrepancies at 90°. This suggests that the free surface condition used in our method can calculate the Rayleigh waves correctly for an epicentral distance less than 90° for a 600 km deep source in this model.

Since a distribution source in space is used in the model (see Section 3.5), the shallowest source depth we can set in the numerical model described in the beginning of this section is 85 km. In order to compare the surface waves excited from a shallower source, we calculated the synthetic seismograms for a 170 km deep source. All other parameters and the numerical model for the calculation are the same as used in the previous calculation. Fig. 6(b) shows the complete seismograms between 10° and 90° in which the PSM and DSM results are superimposed. We see again that the Rayleigh waves in the lateral component (u_θ) are nearly identical for all epicentral distances. However, the amplitude of Rayleigh waves in the radial component (u_r) shows a large difference between the two methods. This means that the accuracy of the surface waves calculated for a shallower source is degraded using the present model. In order to calculate accurately the surface waves using this method, it is necessary to reduce the grid spacing at the free surface in the model (e.g. Rodrigues & Mora 1993).

Much shallower sources, e.g. 20 or 30 km deep, can also be implemented in the model by reducing the grid spacing, but the computation time for such shallow source models will be too long to be carried out on a desktop workstation. The major advantages of this method are that the centre and the inner core can be included in the whole earth model, so that the whole wavefield can be calculated. Since we mainly focus on the core phases that propagate through the inner and outer core of the earth, a 600 km deep source that allows the P and SV waves to be directly radiated into the mantle with a clear separation from the depth phases reflected at the free surface is used in this paper. For the same reason, a deep source (600 km) was also used by Furumura et al. (1998). Because for such a deep source even surface waves can be correctly calculated using the present model (Fig. 6a), we will use the whole wavefield and the complete seismograms obtained by this method to show the effects of arbitrary heterogeneities on seismic wave propagation in the following sections.

The above comparisons for core phases, surface waves and the complete seismograms suggest that it is possible to simulate the wavefield in a 3 D whole spherical earth model with reasonable accuracy by using the method we developed for a 2 D cylindrical earth model. By using the 2 D method, we can significantly reduce the required computer memory and computation time for the whole earth modelling, especially when we wish to simulate seismic wave propagation in the whole earth with arbitrary heterogeneous structures that can be approximated as azimuthally symmetric. Since singularities arise at pole axes θ=0 and π when we solve the equations of motion in 3 D spherical coordinates, it might be difficult to calculate the wavefield around these positions. In 2 D cylindrical coordinates, by avoiding the singularity at the centre of the earth, we can efficiently calculate the whole wavefield for a 360° earth, and obtain a clear image of the wave propagation in the whole earth, which is very helpful in understanding the full process of the whole wavefield evolution.

5 Wavefield and synthetic seismograms for the IASP91 model

In this section, we apply the method to the IASP91 model (Kennett & Engdahl 1991) to simulate wave propagation and the motion at the surface for a 2 D cylindrical radially heterogeneous earth. The results are presented as sequential snapshots of the wavefield and synthetic seismograms along the surface.

The numerical model used in this simulation consists of three subdomains, and the grid numbers in the radial and lateral directions in each subdomain are the same as in the model used in the last section. The source is a double couple line source located at 600 km depth with moment tensor components M_rr=−1.0 and M_θθ=1.0. The time history of the source is a Herrmann function with a pulse width of 30 s. The number of gridpoints per minimum S wavelength for the maximum grid spacing in this model is 2.58, and Δa^ref is 14.0 km, which allows a time interval Δt=0.25 s. The total number of time steps calculated is 12 000 for synthetic seismograms of 3000 s duration.

The anelastic attenuation of the medium is incorporated in the calculation by using the following attenuation coefficients given by Graves (1996):

where Q is the anelastic attenuation factor for a reference frequency f₀. In the calculations, these coefficients are multiplied by the values of stress and velocity at each gridpoint at each time step (Furumura et al. 1998). The Q used here is then strongly frequency dependent, although a flat Q may be preferred in global wavefield modelling. We use this scheme in this calculation because it requires little computer memory and computation time. The Q_P/Q_S effect can also be separately included in the calculation by following the recent studies of Hestholm (1999), Hestholm & Ruud (2000) and Olsen et al. (2000). Efficient implementation of such new techniques will be the subject of a future study.

5.1 Wavefield snapshots

In Fig. 7, the sequential wavefield snapshots at 16 time steps are displayed in order to show the generation and propagation of various seismic phases in the earth. The P and SV waves are represented by red and green colours, respectively, and the contributions from P and SV waves are calculated from the divergence and curl of the wavefield:

The sequential wavefield snapshots give us a clear image of the history of seismic wave propagation inside the whole earth radiated from a double couple source. They are very helpful tools in identifying the generation and evolution of various phases and showing clearly the relative amplitude of each phase at various epicentral distances.

In the frames at 300 and 450 s, artefacts can be seen prior to the arrival of PKI at the centre of the earth. This is caused by the accumulation of very weak noise on gridpoints at the centre. As shown in Fig. 1, the 256 lateral gridpoints are almost located at the same position and the background noise is enhanced more than 200 times, which is clearly visible at the centre. This artefact can be seen in all the following snapshots, which could be reduced by reducing the number of gridpoints around the centre in the display.

At 600 s, the core phase PK propagates in the outer core and PKI travels through most of the inner core, while pPK and SK reach the inner core. P, pP, sP, S, pS and sS propagate in the mantle with clearly separated wavefronts, and P_diff and pP_diff diffracted at the CMB propagate in the lower mantle. In the upper mantle, the surface multiple PP travels with strong energy. At 750 s, PKIK transmits through the inner core into the outer core and pPKI propagates in the inner core. pS and sS penetrate into the outer core and convert to the core phases pSK and sSK. In the mantle, ScS approaches the surface, the sScS phase is generated at the CMB, and the conversion SKS appears near the CMB. In the uppermost mantle, both P and S surface multiples travel with strong energy. At 900 s, PKP and PKIKP propagate out of the core and into the mantle, and PKS is generated at the CMB when PKP transmits from the core into the mantle. pPK and pPKIK propagate in the outer core. The SKK phase, the internal reflection of SK at the CMB, propagates along the CMB in the uppermost outer core. In the mantle, P_diff arrives at the surface, and the P and S surface multiples and reverberations within the upper mantle cause a very complex wavefield pattern. SKS and S_diff are propagating within the lower mantle, and the weak CMB reflection sScS is also visible.

In the last three frames at 600, 750 and 900 s, the PKJ and pPKJ phases with very weak shear wave energy converted from PK and pPK at the inner–outer core boundary can be seen clearly in the inner core. Since the wavefront of PKJ is connected to that of PKI at the inner–outer core boundary and the connection points of the wavefronts run faster than the wavefront of PKJ inside the inner core, the wavefront of PKJ inside the inner core gradually grows to a shape with a reversed convex to that of PKI, as displayed in the frame at 600 s. The wavefront of PKJ has a similar shape to that of PKI only for a short time after PKI propagates into the inner core. The PKJ phase then advances in the inner core with the wavefront convex upwards as shown in the frame at 750 s. In the lower half of the inner core, the upward propagating PKIJ reflected at the inner–outer core boundary has a similar shape of wavefront to PKI. At 750 s, PKIJ is overlapped by the wavefront of a strong pPKI but can be identified with care, while the wavefronts of PKJ and PKIJ form a ‘circle’ shaped wavefront in the lower half of the inner core. The situation of pPKJ is similar to that of PKJ, which is clearly visible in frames at 750 and 900 s. The ‘circle’ shaped wavefront formed by downward propagating pPKJ and upward propagating pPKIJ is very clear in the frame at 900 s. Until now, there have been no confirmed reports of the identification of the PKJKP phase, which propagates through the mantle and outer core as a compressional wave but traverses the inner core as a shear wave. The difficulty in detecting PKJKP is caused by the inefficiency of P to S and S to P conversions at the inner–outer core boundary. In a very recent paper, Deuss et al. (2000) reported the identification of pPKJKP and SKJKP from two deep earthquakes (the 1996 June 17 Flores Sea event at 584 km depth and the 1994 June 9 Bolivia event at 647 km depth) in long period observations (20–30 s). The wavefield snapshots of the whole earth presented here provide us with a clear image to ‘see’ the generation and propagation of these shear waves in the inner core.

5.2 The synthetic seismograms

The synthetic seismograms at the earth's surface are displayed on a reduced timescale extending to 1800 s in Fig. 8. In the lateral component, we see the S phase with large amplitude, the clear core reflection ScS, the SKS and the triplication of these phases at around 80°. Surface multiple (ScS)₂ can be clearly seen extending to between 120° and 140° between SS and SSS. Other core reflections such as PcS and sScS can also be clearly discerned. S and its surface multiples SS and SSS extending to large range with high amplitude are prominent features in the lateral component seismograms. In the radial component, the first major phase is P with large amplitude up to about 100°, followed by a superposition of pP and PP with large amplitude between 30° and 50°. SP is clear after S and S_diff extending in a wide range between 40° and 130°. Core diffracted P_diff is visible beyond 100° up to 120°. Various branches of core phases PKP and pPKP can be seen beyond 120°, and P′P′ is visible between 65° and 85°.

Fig. 9 shows the core phases in the radial component excited from a 600 km deep explosion source for the same model as in Fig. 8 together with the ray theoretical arrival times of these phases for the IASP91 model. Since the whole core is included in the model, we can obtain the complete seismograms for various branches of core phases between epicentral distance of 116° and 244°. Fairly good agreements in arrival times can be seen for all the core phases compared with the ray theoretical traveltimes. Around 144° and 216°, the triplications of the PKP branches have large amplitudes. The triplications with large amplitudes can also be seen for pPKP branches around 148° and 212°. The core diffracted branches AB of PKP and pPKP, which propagate down into the core rather than up into the mantle, extend well over their ray theoretical limits of 170° and 179°, respectively. This was observed by Shearer (1991) for global body wave phases by stacking the long period GDSN seismograms.

The diffracted PKP phase can be seen in long period seismograms and arrives before the PKP(DF) branch. In Fig. 10, we enhance the amplitude of the seismograms in Fig. 9 to show the phases around the triplication at 144° from 1000 to 1250 s and between 125° and 180°. The PKP_diff can be clearly seen ahead of PKP(DF) between 130° and 144°. In the acoustic simulation of P wave propagation in an axisymmetric whole earth model by Thomas et al. (2000), the PKP_diff phase with a dominant period of 11 s is calculated for a source at a depth of 192 km. Comparing our results with their results, we see that PKP_diff is more pronounced in this simulation with a source time function of 30 s width. The waveforms of PKP_diff and the following PKP(DF) overlap in these long period seismograms but are clearly separated in the short period results of Thomas et al. (2000).

We also calculated the complete synthetic seismograms for a shallower source located at a depth of 170 km using the same model as in Fig. 8. The synthetic seismograms are displayed in Fig. 11, where both the time and the amplitude scales are the same as those used in Fig. 8. Comparing Fig. 11 with Fig. 8, we find that the seismograms for the shallow source are dominated by Rayleigh waves and surface multiples that propagate along a shallow zone below the surface. The arrival times for core phases PKP and PKS are very close to those for pPKP and pPKS, respectively, so it is difficult to have a clearly separated waveforms for these core phases. Since this paper mainly focuses on the complete core phases in the whole wavefield, we prefer to use the results from the deep source (600 km) to show the whole wavefield evolution as in the snapshots in Fig. 7 and the effects of arbitrary heterogeneities as shown in the examples in Section 6.

6 Examples for models with heterogeneity above the CMB

In recent years, detailed studies of broad band waveforms for various secondary phases reflected, refracted or diffracted from the CMB have revealed the structures of velocity perturbations above the CMB and its topography (e.g. Lay et al. 1998). In this section, we apply our method to models including a velocity perturbation zone just above the CMB to examine the wavefield anomalies caused by the heterogeneous structure. Although recent studies of CMB structure were usually performed by using short period seismograms to locate the relatively short wavelength heterogeneity (e.g. Ritsema et al. 1998; Vidale & Hedlin 1998), we carry out the modelling for a large scale heterogeneity for long period body waves in order to obtain a clear picture of the wavefield variations inside the perturbation zone. The mechanism of the double couple line source is the same as used in the previous modelling and the source time history is a Herrmann function with a pulse width of 30 s. The source depth is 600 km. The numerical model is also the same as used in the previous modelling. The same anelastic attenuation coefficients as used in the last section are employed in these examples. The ellipse shaped velocity perturbation zone extends 2400 km along the CMB and up to 600 km into the mantle with a 5 per cent reduction in both Pet al. and S wave velocities.

6.1 Example 1

In this example, the centre of the velocity perturbation zone is located 60° from the epicentre. In Fig. 12, the wavefields inside and around the low velocity heterogeneity are displayed at four time steps for both models—including the velocity perturbation and without the perturbation. At 600 s, P_diff can be seen propagating more slowly when the low velocity zone is introduced in the model than in the laterally homogeneous model, pP, sP are travelling inside the heterogeneity and S and ScS have just entered this zone. At 750 s, the sP_diff is seen propagating slowly with larger amplitude in the laterally heterogeneous model. The SKS phase appears inside the low velocity zone, and both SKS and ScS there have larger amplitude for the laterally heterogeneous model, and the pS phase approaches the low velocity zone. At 900 s, it can be clearly seen that SKS and S_diff are propagating more slowly, with larger amplitude in the laterally heterogeneous model. The pS phase is propagating inside the low velocity zone without visible change and sS and sScS approach this zone. At 1050 s, we see that pS_diff propagates slowly with larger amplitude in the laterally heterogeneous model and sS_diff is travelling in the low velocity zone more slowly with larger amplitude than in the laterally homogeneous model. These snapshots give us a clear image of the effects of the low velocity zone above the CMB on the features of various phases propagating through it. In these snapshots, ‘trapped’S wave energy just above the CMB can be seen, which is an artefact associated with the fluid–solid boundary condition. This is the ringing effect due to the Fourier transform at the discontinuity at the CMB. This artefact decays very quickly in the radial direction and does not propagate into the mantle and the outer core, so it does not affect the core phases, as shown in Section 4.2.

The differential seismograms obtained by subtracting the results for the model with velocity perturbation from those without perturbation are shown in Fig. 13. In both the lateral and the radial components, the anomalies of S, ScS and SKS phases are prominent around the triplication at around 80°. Delayed arrival times can be seen for ScS between 60° and 80°, for S over 80° and for SKS between 80° and 100°. The strong anomalies can also be seen within the range where the depth phases sS and pS and the core phases sSKS and pSKS have little arrival time difference. In the radial component, the anomalies for P and P_diff can be seen between 75° and 120° and for pP and pP_diff between 80° and 125°. The differential seismograms show the complete wavefield anomalies caused by the low velocity zone above the CMB. By studying these features, we may perform a detailed analysis of waveforms at certain epicentral distances in order to look for quantitative relations between velocity perturbations and anomalies in arrival times and amplitudes. The differential traveltimes between S, ScS and SKS have been used for locating the spatial distribution of velocity perturbations in the lower mantle (e.g. Ritsema et al. 1998).

6.2 Example 2

In the next example, we put the centre of the low velocity zone at an epicentral distance of 130° to investigate the effects of the velocity perturbation on the core phases. Precursors to the PKP phase in high frequency seismograms (up to 2.5 Hz) have been used for identifying a low velocity region above the CMB (e.g. Vidale & Hedlin 1998). However, this modelling provides us with the wavefield anomalies caused by the low velocity region in the long period range. The differential seismograms in the radial component are displayed in Fig. 14. The delay of arrival times for various branches of PKP and pPKP can be seen clearly between 140° and 150°. The largest delay of arrival time occurs between 144° and 150° for both PKP and pPKP. The arrival time delay for PKP(AB) and pPKP(AB) is visible for epicentral distances greater than 180°. The PKP_diff close to 144° ahead of PKP(DF) is seen to be affected strongly by the heterogeneous zone, and the anomaly of this diffracted phase is visible between 135° and 144°.

7 Discussion and Conclusions

We have developed a method for modelling the seismic wave propagation in a whole earth model by using the Fourier pseudospectral method. With this method, we can simulate the complete wavefield in a 2 D cylindrical whole earth model with arbitrary heterogeneity. The examples for the IASP91 model and models with a velocity perturbation above the CMB show that the features of wave propagation and the effects of heterogeneity on various seismic phases can be clearly seen with sequential snapshots and differential seismograms.

The model used in this method is defined in 2 D cylindrical coordinates rather than as a 3 D spherical earth model. However, the waveforms excited in the 2 D cylindrical model from a ‘line’ source can be mapped to the corresponding waveforms in a 3 D model from a ‘point’ source. The comparison simulations carried out in the previous sections show very good agreement between the 3 D results and the mapped 2 D results, except for a slight difference between a few phases caused by the large grid spacing. The agreement between them would be improved if we used finer grids in the model. For surface waves from a shallow source as shown in Fig. 6(b), it is necessary to use a much finer grid spacing below the free surface to ensure accuracy. The comparison suggests that it is possible to simulate the wavefield in a 3 D spherical earth model by using a 2 D cylindrical model. This will significantly reduce the computer memory and computation time required when we perform forward modelling for whole earth models to understand the full process of wavefield evolution, especially when short wavelength arbitrary heterogeneity is included in the model. The merit of the 2 D cylindrical modelling is that we can calculate the whole wavefield for a 360° earth with arbitrary heterogeneity in order to study the effects of the heterogeneity on the whole wavefield, although the heterogeneity is limited to the great circle cross section. However, the scattering and focusing/defocusing effects from out of plane wavefields cannot be accurately predicted in a 2 D cylindrical model. This method is useful for models in which the heterogeneity can be approximated as azimuthally symmetric.

The problems related to the centre of the earth are also present in the 3 D spherical whole earth model, and the ‘extension scheme’ we proposed to avoid the singularity at the centre, the multidomain technique and the ‘smoothing scheme’ to increase the time interval are also useful for a 3 D spherical model including the centre. The schemes we presented in this method can be directly extended to the 3 D case in order to solve the problems related to the centre of the earth.

Recent studies of the earth's structure have been performed in a high frequency range (period less than 20 s) (e.g. Hedlin et al. 1997; Vidale & Hedlin 1998; Ritsema et al. 1998) to investigate short wavelength heterogeneity. Therefore, it is desirable to simulate the effects of small scale heterogeneity in the earth on the short period global wavefield. In this paper, we calculated the wavefield for a dominant period of 30 s. The simulation required 76 Mbyte memory in a single precision calculation with a CPU time of 7.8 hr on a DEC Alpha workstation ‘au500’ (500 MHz clock speed) for a 3000 s duration of wave propagation. The method can be used for modelling in the high frequency range, but the required computation time will be much longer. For a dominant period of 7.5 s, the CPU time for calculation will be about 1300 hr using the same workstation, which is unreasonably long. Therefore, further improvement to reduce the computation time is necessary to apply the method to the high frequency range.

Acknowledgments

We are very grateful to N. Takeuchi at the Earthquake Research Institute, the University of Tokyo, for use of his DSM program package. Constructive comments by P. Cummins, M. Sambridge and an anonymous reviewer are appreciated. This study was partially supported by the Superplume Project funded by the Science and Technology Agency of Japan.

References

Author notes