Accurate numerical methods for solving the elastic equation of motion for arbitrary source locations

Summary

In our previous studies we derived a general criterion for optimally accurate numerical operators for computing synthetic seismograms. In this paper we extend our results by deriving source representations that are ‘tuned’ to match the errors of the numerical operators, so that optimally accurate synthetic seismograms are obtained. The same accuracy can be obtained whether the sources are at nodes of the numerical grid or between nodes. The accuracy of the synthetic seismograms is confirmed by numerical experiments. These results should facilitate inversion of seismic waveform data for 3-D Earth structure and earthquake source parameters.

earthquake source representation, synthetic seismograms, waveform inversion

1 Introduction

Accurate and efficient methods for computing synthetic seismograms in laterally and vertically heterogeneous media are required to perform inversion of seismic waveform data for 3-D Earth structure (e.g. Geller & Hara 1993). Geller & Takeuchi (1995, 1998, hereafter referred to as GT95 and GT98, respectively) and Takeuchi & Geller (2000, hereafter referred to as TG00) derived optimally accurate numerical operators for computing synthetic seismograms that are well suited to waveform inversion studies. Schemes using these numerical operators require far less CPU time to compute synthetics of any given accuracy than conventional numerical schemes (e.g. an improvement by a factor of about 50 is obtained using optimally accurate finite-difference operators for a 2-D heterogeneous medium; see TG00 for details).

In deriving the above numerical schemes, GT95, GT98 and TG00 considered only the discretization errors of the numerical operators. In this paper we extend their treatment by also considering the discretization error due to the force term. We first derive such source representations for the Direct Solution Method of Geller & Ohminato (1994, hereafter referred to as GO94); however, we also demonstrate that our approach can be applied to the staggered grid and pseudo-spectral methods. The numerical examples presented in this paper confirm that when the theoretical results derived here are implemented using optimally accurate numerical operators, comparable accuracy is achieved for sources in the interval between nodes as for sources at nodes, while the accuracy of synthetic seismograms obtained by conventional source representations seriously degrades for sources between nodes.

2 General Condition for Optimal Accuracy

2.1 Theoretical Background

The weak form of the equation of motion (see GO94) is 

(1)

where ω is the frequency, ρ is the density, C_ijkl are the elastic constants (including the effects of anelastic attenuation), u_i is the displacement, φ^(m)_iis the i -component of the m th trial function, f_i is the external force, V is the entire volume, the superscript * denotes the complex conjugate, and summation over repeated indices is implied. Note that all of the results in this section are in the frequency domain.

Eq. (1)is transformed to a system of simultaneous linear equation by representing the displacement as a superposition of trial functions: 

(2)

The trial functions can be either low-order polynomials (e.g. linear splines), or analytic functions (e.g. those used by Geller & Hatori 1995, hereafter cited as GH95). Trial functions are usually not explicitly defined for weak-form finite-difference methods, but they are implicitly low-order polynomials such as linear splines. In this paper the trial functions are assumed to be defined such that each trial function is zero everywhere except within one grid interval of one particular node. (Note that the grid spacing is not required to be uniform.)

We substitute eq. (2)into eq. (1)to obtain 

(3)

where Tand Hare the mass and stiffness matrices, respectively, dis the vector of expansion coefficients, which are the unknowns in eq. (3)and gis the discretized external force vector (see eqs 38and 39in GO94). The elements of the matrices on the left-hand side of eq. (3)and the force vector on the right-hand side of eq. (3)are as follows: 

(4)

 

(5)

 

(6)

When non-analytic trial functions such as linear splines are used in eq. (2), the numerical operators in eq. (3)will not in general satisfy the general criterion for optimally accurate numerical operators derived by GT95 (). However, as shown by GT95, GT98 and TG00, non-optimally accurate operators can be replaced by the corresponding optimally accurate operators in a straightforward fashion.

The choice between optimally accurate operators of different types depends on the nature of the problem (Mizutani et al. 2000). For example, for homogeneous media pseudo-spectral spatial operators and fourth-order Runge–Kutta temporal operators outperform optimally accurate second-order (in both space and time) finite-difference operators, but the latter are more cost effective for a heterogeneous medium with sharp boundaries between layers. As the Earth is a heterogeneous medium with sharp internal boundaries, optimally accurate local schemes (finite-difference or finite-element) appear best suited to actual seismological applications. The discussion in this paper is therefore limited to such numerical schemes.

2.2 Extension Of Gt95 Theory To The Source Term

GT95 derived a general criterion for optimally accurate numerical schemes (their eq. 2.20). However, their theoretical analysis did not consider the source term. Here we extend GT95's theory by deriving a general criterion that the source term should satisfy to obtain an optimally accurate solution.

The following discussion is an extension of of GT95. We use first-order perturbation theory to estimate the error due to the source term. The exact solution d⁽⁰⁾and the numerical solution dare, respectively, expressed by a normal-mode summation: 

(7)

 

(8)

where d_kis the k th eigenfunction of the numerical matrix operators Tand H, and the superscripts ⁽⁰⁾denote the corresponding exact quantities. We assume that d⁽⁰⁾_kare normalized so as to satisfy 

(9)

where δ_{k ′k}is a Kronecker delta, and that d_kare expressed as follows: 

(10)

Comparing eqs (7)and (8)and using the approximation 1 /(1 +x) ≈ 1 −x, we obtain 

(11)

where δd, δT, δH, δgand δd_kdenote the difference between the numerical and exact quantities.

We now evaluate the n th mode component of the error. In other words, we write eq. (11)in the form 

(12)

and we evaluate δd_n. Using eqs (9)and (10), we obtain 

(13)

where 

(14)

and 

(15)

Substituting the well-known relation (see e.g. Geller et al. 1990) 

(16)

into eq. (13)we obtain 

(17)

Eq. (17)is a general result, but, as discussed by GT95, the contribution of the n th mode to the exact solution and numerical solution will be dominant when the frequency ω is close to the eigenfrequency of the n th mode ω_n. For this case, we can ignore the last term of eq. (17)because its denominator does not approach zero as ω→ω_n. When ω is around the n th eigenfrequency ω_n, we can evaluate the relative error as follows: 

(18)

Note that we use the well-known relation between modal expansion coefficients d⁽⁰⁾_nand excitation coefficients g⁽⁰⁾_n(eq. 2.14 of GT95) in deriving the last term of eq. (18).

In the terms on the right-hand side of eqs (18), the second term will in general become very large when ω→ω_n, i.e. the relative error will degrade significantly in the neighbourhood where we most desire accuracy. To prevent this, we require the following condition to be satisfied: 

(19)

This is the general criterion for optimally accurate operators derived by GT95 (their eq. 2.20).

Let us further consider eq. (18)for the case of optimally accurate operators (i.e. operators that satisfy eq. 19). Substituting eq. (19)into eq. (18), we obtain 

(20)

In the transformation from the second line to the last line of eq. (20), we used eq. (19). From eq. (20), we can see that if the condition 

(21)

is satisfied, we can obtain an optimally accurate solution when optimally accurate operators are used. We assume that the second term on the right-hand side of eq. (21)can be neglected (it is not clear that this will always be the case in practice; this is a topic for further research), in which case eq. (21)reduces to 

(22)

This is the condition for an optimally accurate solution.

We now consider the intuitive meaning of the above derivation. If the numerical operators are not optimally accurate (i.e. do not satisfy eq. 19), then the second term on the right-hand side of the first line of eq. (20)will dominate the error of the numerical solution, and the error due to the force term will be comparatively negligible. However, if we are using optimally accurate numerical operators, then the contribution of the first term on the right-hand side of the first line of eq. (20)is no longer negligible. Eq. (22)shows that in this case the error of the synthetic seismograms can be minimized if the error of the force term is ‘tuned’ to match the error of the numerical operators.

3 Source Representation

We now derive source representations that satisfy eq. (22). We first derive the source representation for a point force at an arbitrary location when we use analytic trial functions. We then apply this representation to the case of numerical trial functions and derive a source representation satisfying eq. (22). Note that the source representations derived in this section are applicable to sources at nodes as well as to sources between nodes. The following discussion is for sources between nodes, but the same formulation can be applied to sources at a node because their discretized representation is the same as that for a source between nodes, which is infinitesimally close to the node. (See the discussion following eq. 48, below.)

3.1 Source Representation For Analytic Trial Functions

As the trial functions used in weak-form solutions are continuous, it is not possible to account for the displacement or traction discontinuity for a seismic source between gridpoints solely in terms of a trial function expansion. Such discontinuities must therefore be accounted for by representing the displacement in the vicinity of the source as the sum of a particular solution and a trial function expansion.

We begin by considering the weak-form results derived by GH95, who used analytic trial functions that are solutions of the homogeneous (source-free) equation of motion. GH95 represent the displacement in the source region (see Fig. 1a), V_S, as the sum of a homogeneous solution, v_i, and a particular solution, u^(p)_i, while the solution in the remainder of the medium, V_R, is given by a homogeneous solution only: 

(23)

where u^(p)_iis the particular solution satisfying the strong form of the equation of motion in the source region, V_s: 

(24)

which satisfies rigid boundary conditions: 

(25)

on the boundary of the source region, S_s.

Figure 1

Source region (V_S) and remainder of medium (V_R) for (a) plane-layered medium and (b) general 2-D medium. Note that V_Rextends in all directions to the outer boundary of the medium. The diagram for the 3-D case is not shown, but it is a straightforward extension of Fig. 1(b).

By substituting eq. (23)into eq. (1), the right-hand side of the weak form of the equation of motion is transformed to the following: 

(26)

where n_j is the outward unit normal vector of S_S. To obtain the weak-form solution, we represent the homogeneous solution as a superposition of trial functions: 

(27)

The explicit form of the trial functions for the plane-layered SH problem and the P –SV problem are given by eq. (22)and eqs (38)–(46)of GH95, respectively. Generally speaking, for both analytic and numerical trial functions the index, n, of the trial functions and expansion coefficients in eq. (27)is usually a pointer to a set of indices (see Appendix A of GO94). For example, for the P–SV problems for the plane-layered medium in Fig. 1(a), we have n →{q, l }, where q refers to the node point, z_q, and l to the component of displacement (e.g. vertical or horizontal), which is non-zero at that node for that particular trial function. The values of the trial functions at the various node points for the P–SV problem are thus given by 

(28)

where δ_ilis a Kronecker delta.

We substitute the trial function expansion (eq. 27) into the equation of motion (eq. 26) to obtain the following system of linear equations: 

(29)

where Tand Hare defined in eqs (4)and (5).

As noted above, GH95 use analytic trial functions that are solutions of the homogeneous equation of motion. When such trial functions are used it is possible to integrate the definitions of the matrix elements by parts so that the separate operators on the left-hand side of eq. (29)are replaced by a single operator (see GH95, eq. 11): 

(30)

where the elements of Aare as follows: 

(31)

Note that the surface integrals in eq. (31)are evaluated over all of the surfaces where the tractions of the trial functions have discontinuities. For the medium shown in Fig. 1(a) these are the boundaries between the homogeneous layers. The elements of the force vector hare given by (GH95, eq. 10) 

(32)

3.2 General Source Representation

The particular solution used by GH95 is required to satisfy the rigid (zero-displacement) boundary conditions given in eq. (25). However, it is useful to be able to compute synthetic seismograms for arbitrary particular solutions, as in many cases particular solutions can most easily be computed for an infinite medium (a ‘whole space’) with radiation boundary conditions. We now extend GH95's derivation to the case of a particular solution, U^(p)_i, with arbitrary boundary conditions on the boundary of the source region, S_S, which satisfies the following inhomogeneous equation of motion: 

(33)

We decompose U^(p)_ias follows into the sum of the particular solution satisfying a rigid boundary condition, u^(p)_i, and a remainder, R_i: 

(34)

Substituting eq. (34)into eq. (23), the displacement u_i is expressed as follows: 

(35)

Note that even if we use a particular solution U^(p)_ithat does not satisfy a rigid boundary condition, the relation v_i =u_i still holds at every node.

We now derive some of the properties of R_i. Subtracting eq. (24)from eq. (33)and using eq. (34)we obtain 

(36)

On the other hand, from eqs (25)and (34)we obtain 

(37)

Thus R_i is a solution of the homogeneous equation of motion eq. (36)in the source region V_Rthat satisfies the boundary condition given by eq. (37).

From eqs (36), (32)and (34)we obtain the weak form of the equation of motion: 

(38)

We rewrite eq. (38)in matrix form as follows: 

(39)

where T, Hand Aare defined in eqs (4), (5)and (31), respectively, and the elements of the vector hare given by 

(40)

To evaluate the elements of hwe must know the explicit value of R_k,l. As the trial functions span the space of solutions of the homogeneous equation of motion, R_i can be expressed in terms of a superposition of trial functions. Since each trial function has unit value for only one component at one node and is zero for all other components (see eq. 28), the coefficients of this expansion can be determined essentially by inspection. Using eq. (37)and denoting the discretized vectors corresponding to the values of R_i and U^(p)_iat the nodes by Rand U^(p), respectively, we have 

(41)

We substitute the trial function expansion of R_i into the second term on the right-hand side of eq. (38), and use eq. (31)to obtain the following result: 

(42)

where T′ and H′ are mass and stiffness matrix operators, respectively, defined as in eqs (4)and (5), except that the volume of integration is limited to V_Srather than taken over the entire medium. Similarly A′ is defined as in eq. (31), except that the surfaces of integration are limited to S_S.

Finally, we transform the surface integral in the last line of eq. (42)into a volume integral that is expressed in terms of standard matrix elements: 

(43)

The starting point of the above transformation is the last line of eq. (42). Even if there is a discontinuity in elastic properties across S_S, C_ijklU^(p)_k,ln_j will be continuous. Therefore, rather than evaluating the surface integral on the ‘inner surface’ of S_S, with n_l pointing outward, we evaluate the surface integral on the outer surface, which we designate as S⁺_S, with the normal vector pointing inward. As there are no sources outside V_S, U^(p)_imust be a solution of the homogeneous equation of motion outside V_S, and therefore can be represented as a superposition of trial functions. (If U^(p)_iis not defined outside V_swe can obtain its value outside V_sby analytic continuation, e.g. Carrier et al. 1966, pp. 63–66). We then use eq. (31)to transform the surface integral in the second line of eq. (43)into the volume integrals in the third line. Note that the volume of integration is V_R; i.e. the source region V_Sis excluded. Finally, in the last line of eq. (43)we express the final result in terms of matrix elements. The matrix elements T″ and H″ are defined following eqs (4)and (5). However, the elements of the matrices T″ and H″ are set to zero except for elements in rows corresponding to trial functions for which one of the components has unit value on S_S. Note that the boundary conditions for U^(p)_ioutside the region for which the elements of the matrices T″ and H″ are non-zero are arbitrary.

3.3 Source Representation For Numerical Trial Functions

Now we apply the above results to the case of numerical trial functions. When using numerical trial functions such as linear splines, some of the terms in the above derivation, such as the surface integrals in the first two lines of eq. (43), may not, strictly speaking, be computable. Nevertheless, the above derivation may be regarded as being valid in the sense of the theory of distributions, as the trial functions could be represented as the limit of a series of functions for which the required computations are rigorously possible. Therefore, the last line of eq. (43)may be used in numerical computations even for the case of trial functions for which the first two lines of eq. (43)cannot be rigorously evaluated. See, for example, Strang & Fix (1973, pp. 11–16) or GO94 () for further discussions of this point.

However, we must obtain T″ and H″ in eq. (43)so that eq. (22)is approximately satisfied. (Note that for this case g⁽⁰⁾_nand δg_n in eq. (22)should be replaced by h⁽⁰⁾_nand δh_n, respectively.) Here, we show that eq. (22)is approximately satisfied if we use optimally accurate operators for T″ and H″ on the right-hand side of eq. (43). Expressing eq. (43)in terms of the normal-mode basis (not the modes of the source region V_sbut rather those of the whole region V) and evaluating the error, we obtain 

(44)

On the other hand, because T, H, T″ and H″ are all optimally accurate operators, we can expect that 

(45)

Eq. (45)would be almost exactly satisfied if the number of elements per wavelength (gridpoints per wavelength) is approximately constant over the whole volume. Using eqs (44)and (45), we see that eq. (22)is satisfied. The rigorous quantification of the error when the grid spacing (expressed in terms of elements per wavelength) is not constant remains a topic for future research, but the numerical tests presented in this paper suggest that satisfactory accuracy is achieved for practical applications.

3.4 Example

As an example, we present the explicit elements of the source vector defined in eq. (43)for the case of a homogeneous 1-D string (see GT95) for which the density is ρ, the elastic constant is κ and with a constant grid interval, Δz. The source is located between the n th and (n + 1)th gridpoints. (Note that in this subsection n denotes a particular node rather than a trial function.)

The explicit forms of the matrices T″ and H″ in eq. (43)are as follows when optimally accurate operators are used: 

(46)

 

(47)

The only non-zero elements in T″ and H″ are in the n th and (n + 1)th rows, because only the trial functions corresponding to these rows have non-zero values on S_s(at the nodes z =n Δz and z = (n + 1) Δz, respectively). Thus U^(p)must be evaluated at z = (n − 1) Δz, z =n Δz, z = (n + 1) Δz and z = (n + 2)Δz, because there are non-zero elements in the columns corresponding to these nodes. Using the operators of eqs (46)and (47), we evaluate the non-zero elements of the source vector (eq. 43): 

(48)

where U^(p)_α=U^(p)(αΔz).

For completeness, we discuss the case of sources at nodes, as well as those between nodes. Suppose U^(p)_nis a particular solution for the point force f =δ(z −z₀), where z₀=n Δz +ε is located an infinitesimal distance from the n th node. By expanding eq. (48)in a Taylor series and using the equation of motion satisfied by the particular solution, it can be readily shown that in the limit as ε→ 0 the same solution (to O (Δz²)) is obtained as for the case of three point force sources at the nodes z = (n − 1)Δz, z =n Δz and z = (n + 1)Δz with weights of 1/12, 10/12 and 1/12, respectively. Similarly, for a point force and z₀=n Δz +ε (ε→ 0), we can readily show that the solution obtained using eq. (48)is the same (to O (Δz²)) as the solution obtained for dislocations (displacement discontinuities) with weights of 1/12, 10/12 and 1/12 at z = (n − 1)Δz, z =n Δz and z = (n + 1)Δz, respectively. It can also be shown that the source representation of Cummins et al. (1994) gives the same solution to O (Δz²) when optimally accurate operators are used.

When the source is in a boundary element (i.e. an interval for which one of the nodes is at the outer boundary), T″ and H″ are defined in a similar fashion to the boundary elements of the matrix operators Tand H. For example, when the source is between the zeroth and first gridpoints, T″ and H″ are defined as follows: 

(49)

 

(50)

Thus, for this case the non-zero elements of the source vector (eq. 43) are defined as follows: 

(51)

 

(52)

As the computational domain required for the evaluation of the non-zero elements of his small compared with the whole region, the CPU time required to evaluate his negligible.

4 Time Domain Solutions

In the above section, we derived results in the frequency domain. We now derive the corresponding results for the time domain. Note that the time domain solution presented here is for the weak form of the equation of motion (e.g. the weak-form finite-difference method, GT98 and TG00). Such methods differ from staggered grid finite-difference methods (e.g. Virieux 1986) in the following two points. (1) The dependent variable is displacement (rather than stress and velocity). (2) Boundary conditions (such as free surface conditions and continuity of traction) are natural boundary (continuity) conditions that are automatically satisfied (see GO94), whereas they must be imposed explicitly for the staggered grid methods.

4.1 Equation Of Motion For Arbitrary Sources

The equation of motion in the time domain is as follows: 

(53)

where the temporal derivative operator, X, the spatial derivative operator, K, and the external force term, p, are, respectively, the time domain counterparts of −ω²T, −Hand hin eq. (29). Note that GT98 and TG00 used Ato denote the temporal derivative operators; here, however, we use Xto avoid confusion with Ain eq. (30). The explicit forms of the optimally accurate operators for Xand Kare given by GT98 and TG00, and in this section cis the discretized displacement: 

(54)

where m is a pointer to the set of subscripts for (m_x, m_y, m_z), Δx, Δy, Δz are the spatial grid intervals, and Δt is the temporal grid interval. Hereafter we write u^N_mrather than u (m_x Δx, m_y Δy, m_z Δz, N Δt).

The source vector p, which is the time domain counterpart of eq. (43), is given by the following: 

(55)

where X″ and K″ are the temporal and spatial derivative operators for the rows for which the corresponding node is inside the source region V_sor on its boundary. The definition of the source region is the same as for the frequency domain case (see Fig. 1). U^(p)is the discretized (in both space and time) particular solution. Note that the boundary conditions for the particular solution are arbitrary. The discretization of the particular solution is the same as in eq. (54).

4.2 Example

As an example, we consider a homogeneous 1-D medium with density ρ, elastic constant κ and spatial and temporal grid intervals Δz and Δt, respectively. We assume the source is located between the m th and (m + 1)th spatial gridpoints. The explicit form of the non-zero elements of pin eq. (55)is as follows when optimally accurate operators (see GT98) are used: 

(56)

Thus the elements of the source vector for the time domain case are non-zero only for the elements that correspond to nodes on S_s. We assume [U^(p)]⁻¹_m= 0, which is required in evaluating p⁰_m. Note that, if we want to evaluate p^N_mto the N_maxth time step, we need to evaluate [U^(p)]^N_mfor each time step up to the (N_max+ 1)th time step. When the source is in a boundary element, X″ and K″ in eq. (55)are defined in a similar fashion to the boundary elements of Xand K, by analogy to eqs (49)and (50). For examples, if the source is between the zeroth and first spatial gridpoints, the non-zero elements of pare defined as follows: 

(57)

4.3 Comparison With The Source Box Approach

The above time domain solution is closely related to the source box approach of Alterman & Karal (1968). Although we do not present a derivation here, Takeuchi (1996) showed for frequency domain numerical solutions that, if we use the analytic solution for an infinite medium as the particular solution and use the conventional numerical operators rather than the optimally accurate numerical operators, the resulting formulation is equivalent to the source box approach. Essentially, the same derivation can be applied to the time domain case as well.

The methods of this paper are preferable to the source box formulation for several reasons. The first is the improvement in accuracy due to using optimally accurate operators (although optimally accurate operators could also be used in the source box formulation). A second reason is the simplicity of coding. The source box approach solves two coupled equations: the equation of motion outside the source region and that inside the source region. These two equations use different dependent variables, and the boundary conditions relating these variables at the boundary of the source region are essential boundary conditions, i.e. boundary conditions that must be imposed explicitly (see GO94). The software must thus be coded to account for different source locations. On the other hand, in our approach there is only one equation of motion (eq. 39in the frequency domain or eq. 53in the time domain). An arbitrary source is accounted for by the source term hin eq. (39)or pin eq. (53). If the source location is changed, all that is necessary is to change this term. This greatly reduces the complexity of the software, as compared with that required for the source box formulation.

5 Extension To Other Methods

5.1 Staggered Grid Solution

Extending our derivation to the staggered grid approach is straightforward. As an example, we present the explicit elements of the source vector for the case of a homogeneous 1-D string. The discretized staggered grid equations of motion (Virieux 1986) are given as follows: 

(58)

where is the velocity at the m th spatial gridpoint and the temporal gridpoint, and is the stress at the spatial gridpoint and the N th temporal gridpoint.

If the source is between the n th and (n + 1)th spatial gridpoints, we have to evaluate the particular solutions v^(p)and σ^(p)between the and spatial gridpoints. The non-zero components of the force term in eq. (58)are as follows: 

(59)

Here again, the only non-zero elements of the force vector are the elements for the m th and (m + 1)th gridpoints. It can be shown that the source box approach for staggered grid formulation (e.g. Virieux 1986) is equivalent to the above method.

Numerical examples (not presented here) show that the synthetic seismograms for the source between nodes computed by using the force term of eq. (59)in the staggered grid formulation give almost equivalent accuracy to that for a source at the nodes. However, as shown by GT98, computationally tractable optimally accurate operators cannot be defined for staggered grid methods. Thus, staggered grid solutions will be much less accurate than solutions using optimally accurate operators for the same CPU time.

5.2 Pseudo-Spectral Method

The derivations for the source representation in this paper are for a purely local basis, but we show below that they can also be applied to the pseudo-spectral method (PSM). The equation of motion for the PSM is as follows: 

(60)

where X^PSM, K^PSMand pare A′, K′ and fof eq. (22)of Mizutani et al. (2000), respectively. The force term corresponding to eq. (55)is 

(61)

where X″^PSMand K″^PSMare defined in a similar fashion to X″ and K″ in eq. (55).

Eqs (55)and (61)are formally equivalent, but the cost for numerical evaluation of these expressions is quite different. As the PSM uses a Fourier basis, evaluation of eq. (61)requires the evaluation of U^(p), X″^PSMand K″^PSMfor all gridpoints. Furthermore, the evaluation of X″^PSMand K″^PSMis not straightforward. Thus, the evaluation of eq. (61)would require considerable CPU time. However, this difficulty can readily be eliminated by using eq. (55)directly instead of eq. (61). In other words, we use the finite-difference method to evaluate the source term, but we use the PSM to solve the equation of motion. In a PSM formulation using both velocity and stress as dependent variables, we can use the source term in the same fashion for the staggered grid finite-difference method (e.g. eq. 59). Nonetheless, as pointed out by Mizutani (2000), the poor cost–performance ratio of the PSM for heterogeneous media with sharp discontinuities is a strong argument against using the PSM for modelling wave propagation in realistic Earth models.

6 Application To Cmt Inversion

One of the prime reasons for developing accurate and efficient methods for computing synthetic seismograms is to use them to conduct accurate and efficient waveform inversion for Earth structure and earthquake source parameters. In this section, we show how one can use our source representation (eq. 32or eq. 43) in CMT inversion. The quantities that must be computed for CMT inversion are the synthetic seismograms u_i and their partial derivatives δu_i with respect to the perturbation of the source parameters (centroid and moment tensor) at the receiver position r.

The efficient CMT inversion algorithm of Hara (1997) is as follows. First we obtain zby solving 

(62)

where ^†denotes the conjugate transpose of a matrix, * denotes the complex conjugate of a vector and 

(63)

The physical meaning of zis a back-propagated wavefield excited by a point force directed in the i th direction acting at the position of a particular station, r. See of Geller & Hara (1993) for details. Then we obtain u_i(r) by computing 

(64)

and we obtain the partial derivatives δu_i(r) by computing 

(65)

The key advantage of this algorithm is that, as zis common between eqs (64)and (65), the required number of solutions of eq. (62), which is the most computationally intensive step, can be reduced. Using this algorithm only one LU-factorization and one forward- and back-substitution per receiver are required (see Hara 1997, for details).

When we apply the source representation derived in this paper to the above algorithm, there will be no difficulty in defining the terms on the right-hand side of eqs (62)and (64). These are a single force at the receiver and the source term of the actual earthquake, respectively, and we can apply eqs (32)and (43)as they are.

The evaluation of δgin eq. (65)will depend on whether we are using eq. (32)or eq. (43). We can compute δhfor hof eq. (32)as follows: 

(66)

where δu^(p)_k,lis the perturbation of the particular solution caused by the perturbation of the source term.

7 Numerical Examples

In this section, we present numerical examples that confirm that optimal accuracy is achieved if and only if optimally accurate numerical operators are used in combination with the source representation derived in this paper.

Fig. 2 shows the accuracy of synthetic seismograms for various source locations. We consider a heterogeneous 1-D medium (a ‘string’) for which the density ρ (g cm⁻³) and velocity v (km s⁻¹) at a depth of r (km) are given, respectively, by 

(67)

 

(68)

where R = 5000 and 0 ≤r ≤R. We consider point dipole sources at various source locations r₀between 3000 ≤r₀≤ 3025. We compute synthetic seismograms using: (1) optimally accurate operators and a conventional source representation (given by eq. 6) and (2)optimally accurate operators and the source representation derived in this paper (given by eq. 43). We use optimally accurate operators for T″ and H″ in eq. (43)except where otherwise stated.

Figure 2

The error of synthetic seismograms computed by the source representation derived in this paper (solid line) and the conventional source representation (broken line) for various source locations. Vertical lines show the location of numerical nodes. The error of synthetics computed for sources at nodes by the source representation of Cummins et al. (1994) (crosses) is also shown. Optimally accurate operators are used for all cases. The vertical axis is on a log scale.

We use linear spline trial functions for the conventional source representation and use 1001 nodes (i.e. a grid interval Δr = 5 km). To evaluate the error, we compute highly accurate synthetic seismograms using 10⁵grid intervals (i.e. Δr = 0.05 km) and treat them as ‘exact’ solutions. We change the location of the point source such that r = 3000.00, 3000.05, 3000.10, …, 3024.95 and 3025.00; thus the source is always at a node in computing the ‘exact’ synthetic seismograms. To express the relative position of the source location r₀and the surrounding nodes r₁and r₂, we introduce the parameter ξ, which is defined as 

(69)

The particular solution U^(p)required to compute the source representation in eq. (43)is evaluated numerically using optimally accurate operators with ten times finer grids than used in the rest of the medium, with an additional node at the source (Δr ≈ 0.5 km). However, we have confirmed that the method used to evaluate U^(p)does not significantly affect the results as long as it is sufficiently accurate. Fig. 2 shows that using our source representation the accuracy is excellent for all values of ξ, whereas using the conventional source representation the accuracy depends highly on ξ and degrades seriously unless ξ is very close to 0.5.

When the source is exactly at a node (i.e. when r₀= 3000, 3005, …, 3025 for this case), we can use the source representation of Cummins et al. (1994), which explicitly imposes the discontinuity condition at the source. The accuracy of synthetic seismograms computed using this source representation is also shown in Fig. 2. This confirms that our source representation achieves comparable accuracy to this method even when the source is not located at a node.

Fig. 3(a) shows a comparison of the accuracy of the various methods as a function of the number of grid intervals. Because the accuracy of the synthetics is highly dependent on ξ for the conventional source representation, we change the source location for each grid interval to satisfy ξ= 0.283 around r₀= 3000 km. The highly accurate synthetic seismograms used to evaluate the accuracy are computed by the same method as that for the numerical examples in Fig. 2 except that we put an additional node at the source location. As expected, synthetics computed using a combination of conventional operators and the conventional source representation have the worst accuracy. However, if we replace the conventional matrix operators by optimally accurate operators but use the conventional source representation, the accuracy of the synthetics does not improve significantly. If and only if we use both optimally accurate matrix operators and the source representation derived in this paper, we obtain high accuracy. As shown by Fig. 3(a), this combination yields O (Δr²) accuracy, whereas the combination of optimally accurate matrix operators and a conventional source representation gives only O (Δr) accuracy for a source between node points. This is because the conventional source representation of eq. (6)has only O (Δr) accuracy. For this computation (iii in Fig. 3a), the source representation in eq. (43)was evaluated using the optimally accurate operators for T″ and H″. In Fig. 3(b), we evaluate the accuracy when we use the conventional operators for T″ and H″. Even for this case the accuracy is proportional to O (Δr²) and is greatly improved compared with cases i and ii in Fig. 3(a), for which the conventional source representation is used. However, as expected from the theoretical considerations in and eqs (44)and (45), better accuracy is obtained by using the optimally accurate operators for T″ and H″ in eq. (43).

Figure 3

(a) (i) The error of synthetic seismograms computed using a combination of the conventional operators and the conventional source representation (dashed line); (ii) error of synthetic seismograms computed by a combination of the optimally accurate operators and the conventional source representation (broken line); and (iii) error of synthetic seismograms computed by a combination of the optimally accurate operators and the source representation derived in this paper (solid line) for various grid intervals. The vertical and horizontal axes are both on log scales. (b) (iii) the same as (iii) in (a); (iv) the same as (iii) in (a), except that conventional operators are used for T″ and H″ to evaluate the source representation in eq. (43).

The above results also hold for more realistic earth models. Fig. 4 shows synthetic seismograms computed for various combinations of operators and source representations. All of the synthetic seismograms are computed for a spherically symmetric model using the Direct Solution Method software of Cummins et al. (1994). The synthetics are for the toroidal component of velocity seismograms for the isotropic PREM model (Dziewonski & Anderson 1981) for a vertical dip point source at 597 km depth with a step-function time history. The distance between the epicentre and the station is 108°. Computations are made for the period range 20–5120 s, and a six-pole low-pass filter with a corner of 31.25 mHz is applied. We use 200 nodes in the radial direction in the mantle. The very accurate synthetic seismograms used to assess the accuracy are computed using the optimally accurate operators for very fine grids (3200 nodes) with an additional node at the source location. For the computations in Figs 4(a)–(d), the source is located between nodes. The synthetics in Figs 4(a) and (c) are computed using eq. (3)for the conventional and optimally accurate operators, respectively. Because the trial functions are continuous between nodes the displacement discontinuity required at the source cannot be reproduced by the numerical solution; the relatively large errors in Figs 4(a) and (c) are thus expected. Figs 4(b) and (d) show numerical solutions obtained using the source representations derived in this paper (eqs 29and 43) for the conventional and optimally accurate operators, respectively. These figures show that accurate results are only achieved when the source representation derived in this paper is used in combination with optimally accurate operators.

Figure 4

Comparison of the accuracy of synthetic seismograms computed using various combinations of operators and source representations for PREM. (a) Conventional operators and a conventional source representation, (b) conventional operators and the source representation derived in this paper, (c) optimally accurate operators and a conventional source representation, (d) optimally accurate operators and the source representation derived in this paper, (e) optimally accurate operators and a source representation obtained by putting an additional node at the source location and explicitly imposing the discontinuity boundary condition. In each figure, the top trace is the synthetic seismograms for the respective method, the second trace is an almost exact synthetic seismogram (computed using a very fine grid), and the bottom trace is the residual (the difference between the synthetic seismogram and the almost exact seismogram).

Although the source is not at a node for the examples in Fig. 4, we could also have computed the synthetics by placing an additional node at the source and explicitly imposing the discontinuity condition using the method of Cummins et al. (1994), as shown in Fig. 4(e). The accuracy obtained in Fig. 4(e) is comparable to that in Fig. 4(d). Similar accuracy could be obtained if we varied the numerical grid spacing to coincide with the source location. However, since, as discussed below, there are good reasons not to reconfigure the grid for each seismic source, and since the accuracy in Figs 4(d) and (e) is comparable, the use of the methods of this paper appears to be preferable.

8 Discussion

Geller & Hara (1993) showed that the most computationally intensive step in iterative waveform inversion for earth structure is the LU-factorization and the forward and back substitutions in eq. (29), and formulated an algorithm that minimized the required number of such operations. Using their algorithm, the required number of LU-factorizations is only one for each iteration (for each frequency) if the matrix elements are source-independent, while the required number of LU-factorizations would be equal to the number of events if the matrix elements were source-dependent. Thus the lack of dependence of the matrix elements on the left-hand side of eq. (29)on the source location is a key point for efficiently utilizing the algorithm of Geller & Hara (1993). If it was necessary to redefine the grid for each source, as was done to compute the synthetics in Fig. 4(e), this would no longer be possible.

The fact that the matrix elements in our formulation do not depend on the source location also facilitates efficient CMT inversion using Hara's (1997) algorithm, as only one LU-factorization is required for each frequency if the matrix elements are independent of the source location, while the required number of LU-factorizations would be the product of the number of events and the number of iterations if the matrix elements were dependent on the source location. Hara (1997) showed that his algorithm is also useful for simultaneous inversion for CMT solutions and 3-D Earth structure because the LU-factorization operations are common between the CMT inversion and the 3-D structure inversion. However, this also only holds if the matrix elements are independent of the source location. The weak-form solution for arbitrary source locations derived in this paper will thus facilitate efficient and accurate waveform inversion for earthquake source parameters and 3-D Solid Earth structure.

Acknowledgments

We thank Hiromitsu Mizutani for a critical reading of earlier drafts. Discussions with Hiromitsu Mizutani and with Nobuyasu Hirabayashi helped us to clarify the theory in and eqs (44)and (45). We also thank two anonymous reviewers for helpful comments. This research was partly supported by grants from the Japanese Ministry of Education, Culture, Sports, Science and Technology (nos 12740257 and 14540393) and by the Japanese Solid Earth Simulator Project.